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CHAPTER  I 
INTRODUCTION 


Most  design  techniques  for  sampled-data  or  discrete-time  feedback 
systems  are  borrowed  from  those  developed  for  continuous  systems.  For 
example,  the  w-transform  or  fictitious  frequency  method  Is  simply  an 
extension  of  the  Bode  approach^.  Therefore,  the  compensators 
synthesized  by  this  method  are,  in  general,  the  recursive  type. 

It  is  known  that  the  recursive  type  compensators  have  the  following 
f2  Tl 

Inherent  shortcomings  : (1)  stability  problems;  (2)  phase  charac- 

teristic problems;  and  (3)  arbitrary  approximation  problems.  In  recent 
years,  nonrecursive  type  compensators  have  become  increasingly  more 
important.  By  using  the  nonrecursive  type  filters,  one  can  not  only 
eliminate  the  disadvantages  of  the  recursive  filters  but  also  have  a 
simpler  design  procedure  and  an  easier  realization  process. 

This  report  Is  concerned  with  nonrecursive  compensator  design  from 

fundamental  principles  to  detailed  design  procedures.  The  method 

developed  In  this  report  utilizes  the  Fast  Fourier  Transform  (FFT) 

m 

algorithm  and  Is  minicomputer  oriented.  The  simple  FFT  algorithm 
and  the  availability  of  powerful  minicomputers^  enable  us  to  deal  with 
complicated  systems  quite  easily  and  accurately. 

There  are  three  types  of  design  methods  available  for  the 

design  of  nonrecurslve  compensators: 


(1)  Window  method  - The  main  contributors  are  Kaiser 


t8] 


Helms 


[9.10] 


Gold  and  Rader 


[ill 


Rablner^’^,  etc. 


(2)  Analytic  method  - This  method  entails  the  use  of  analytic 
techniques  analogous  to  the  Butterworth  and  Chebyshev 
methods  for  designing  recursive  filters, 
mainly  responsible  for  this  development. 


[12] 

Hermann*  was 


(3)  Successive  approximation  method  - This  efficient  method 

Is  computer  oriented.  Parks  and  McC1ellan^,3\  Rabiner, 

[6] 

Gold  and  HcGonegal  have  developed  a useful  cataioque 
of  compensator  designs. 

Since  the  window  method  Is  simple,  straightforward  and  yet  affords 
plenty  of  room  for  Improvement,  we  will  concentrate  on  the  window  method. 
The  principle  of  the  window  method  Is  based  on  reducing  Gibbs'  phenom- 
enon^ ^ and  improving  the  convergence  of  the  truncated  Fourier 
series  In  the  frequency  domain  along  with  the  interpretation  of  the 
generalized  Gibbs'  pheonomenon.  However,  In  the  literature,  Gibbs' 
phenomena  in  the  time  domain  are  considered  In  detail  while  only  men- 
tioned In  the  frequency  domain.  For  clarity,  we  will  derive  the  time 
domain  Gibbs'  phenomenon  as  a basis  and  develop  Its  frequency  domain 
counterpart  In  a systematic  manner  In  Chapter  II.  The  latter  is  neces- 
sary to  understand  the  phenomena  In  both  domains  before  firmly  estab- 
lishing window  techniques. 

In  Chapter  III,  the  generalized  Gibbs'  phoenmenon  is  established. 

The  trapezoidal  window  function  Is  used  as  a basic  unit  from  which  a 
family  of  new  windows  are  developed,  and  a thorough  understanding  and 
insight  of  windowing  Is  obtained.  Wlndows^®*^’^ , such  as  Hamming, 
Hanning,  Kaiser  etc.,  are  empirical  In  nature;  however,  a theoretical 
base  will  be  presented  In  this  report.  In  the  frequency  domain,  the 
main  lobe  and  side  lobes  of  a windowed  Fourier  series  can  be  clearly 
seen.  Both  the  derivation  and  Interpretation  presented  here  are  believed 
to  be  new. 

[2  3 20 1 

Chapter  IV  concentrates  on  nonrecurslve  filters  . The  most 

common  complaint  in  the  nonrecurslve  filter  field  Is  that  frequency 
transformations  cannot  be  applied;  since  through  their  application  the 
resulting  filter  becomes  recursive.  We  will  develop  a complete,  novel 
technique  from  which  high-pass,  band-pass,  and  band- reject  I on  filters 
can  be  designed  via  frequency  transformations. 
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Chapter  V concentrates  on  the  analysis  of  discrete-time  control 
systems.  The  procedure  for  using  the  FFT  and  minicomputers  to  obtain 
the  frequency  response  and  transient  response  of  discrete-time  control 
systems  will  be  presented  logically  utilizing  block  diagram  transforma- 
tions. Particular  emphasis  will  be  placed  on  the  theorems  associated 
with  the  z-transform.  For  example,  the  Ragazzini-Zadeh  identity^2^  Is 
repeatedly  applied  to  simplify  and  manipulate  the  system  block  diagram. 

A unique  and  novel  stability  study  based  completely  on  return  difference 
is  utilized  to  derive  a stability  criterion  which  is  believed  to  be  new. 


The  design  of  discrete-time  control  systems  is  investigated  in 
Chapter  VI.  Specifications  of  example  closed  loop  control  systems  are 
presented  and  analyzed  using  the  FFT,  Fourier  Analyzer,  and  window 
techniques.  The  design  procedure  through  which  realizable  nonrecursive 
compensators  are  designed  utilizing  the  window  method  Is  presented 
in  detail.  Two  frequently  used  examples  of  discrete-time  control  sys- 
tems (type  "0"  and  type  "1")  are  discussed  in  detail.  The  resulting 
compensators  are  synthesized  and  their  performance  verified. 


In  summary,  this  report  presents  an  FFT-orlented,  minicomputer- 
aided,  design  procedure  for  nonrecursive  compensator  synthesis  of  discrete- 
time feedback  systems.  The  method  is  simple,  powerful,  and  more  accurate 
than  previously  existing  methods. 


CHAPTER  II 

THE  GIBBS'  PHENOMENON 
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2.0  INTRODUCTION 

It  is  well  known^'^  that  for  a periodic  function  with  bounded 
variation,  the  Fourier  series  converges  at  every  point,  particularly  at 
points  of  discontinuity.  In  1899  Gibbs  pointed  out  that  convergence 
of  the  truncated  Fourier  series  behaves  In  a quite  different  manner 
at  points  of  discontinuity.  Bocher,  in  1906,  extended  Gibbs'  result 
and  established  a theorem  commonly  called  the  Gibbs'  phenomenon. 

In  this  report,  the  phenomenon  originally  discovered  by  Gibbs  and 
BOcher  is  referred  to  as  the  Gibbs'  phenomenon  in  the  time  domain. 
However,  the  same  phenomenon  exists  in  the  frequency  domain  and  will  be 
referred  to  as  the  Gibbs'  phenomenon  in  the  frequency  domain.  Gibbs' 
phenomenon  is  an  important  aspect  associated  with  analysis  and  design 
of  nonrecursive  filters  and  compensators. 

Let  us  start  with  analysis  of  the  truncated  Fourier  series. 

2.1  TRUNCATED  FOURIER  SERIES 

Let  f(t)  be  a real  periodic  function  of  time  t with  period  T and 
bounded  in  (-  j,  j) . The  Fourier  series  expansion  of  f(t)  is  given  by 


f(t) 


-t 


k*»-® 


ake 


Jk^Qt 


(1) 


i 


I 


where 


"j  ktV- 


For  computing  advantages  and  practical  purposes,  wc  often  take  only 
a finite  number  of  terms  of  the  Fourier  series  in  Eq.  (I).  In  other 
words,  we  calculate  the  partial  sum  instead  of  dealing  with  the  entire 
infinite  series.  Consider  the  nth  partial  sum  (which  contains  ?n+l 
terms) : 


■ £ 


jki»0t 


Substituting  Eq.  (2)  into  Eq.  (4)  yields 


d (t-r.)  f(r.)  dr. 
n 


where 


jk«..0t 


Equation  6 is  called  the  periodic  Dirichlet  kernel1  1 or  window  func- 
tion in  the  time  domain.  With  some  algebraic  manipulations,  it  is  easy 
to  arrive  at  the  following: 


dn(t) 


w0  S,n  l("  + 2}  VJ 


2 it  si 


■m 


According  to  Eq.  (5),  the  partial  sum  of  a Fourier  series  Is  given  by 

the  convolution  of  f(t)  and  d (t)  over  one  period.  The  Interpretation 

n 

of  Eq.  (5)  for  the  partial  sum  of  Eq.  (I)  Is  straightforward;  however 
it  is  not  simple,  since  the  window  function  [Eq.  (7)]  is  complicated. 


A second  interpretation  In  terms  of  the  Fourier  Transform  results 
if  we  let  F(Jw)  be  the  Fourier  Transform  of  f(t),  then 


L Hi“ 


) e^wt  du 


From  Eq.  (I),  we  have 


F(jw)  - 2it  ^ ak  5 (w-kwg) 
k«-°° 

The  nth  partial  sum  of  f(t)  Is  given  by  Eq.  (A),  thus,  the  Fourier 

Transform  of  f (t)  Is 
n 


Fn(jw)  ■ 2 ir  a^  6(w“ku>0) 
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In  general,  we  have 


V*> 


57  f Fn(Jw)  d“ 

/-« 


f% 

2*  J-™,. 


F(ju)  c*wt  du 


However,  from  the  definition  of  the  Fourier  Transform, 

F(»  - j 


f(t)  e-Jut  dt 


(ID 


(12) 


and  substituting  Eq.  (12)  Into  Eq.  (11)  using  the  dummy  variable  x for 
t In  Eq.  (12) , yields 


, . , f n“°  , -Jux  Jut 

U)-xLL,M'  * 


dx  du 


fc  /'  '<»>  J ° 

1 J-m  J- nuA 


•x) 


du  dx 


Eq.  (13)  can  be  rewritten  as 


f (t)  - I f(x)  w (t-x)  dx 

n J-m  " 


where  w (t)  Is  called  the  aperiodic  Dlrlchlet  window, 
n 


(13) 


(1M 
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sin  (nu^t) 
irt 


(15) 


(16) 


Now,  we  recognize  that  the  aperiodic  Dlrlchlet  window  Is  an  ideal  low 

pass  filter  with  cutoff  frequency  nu^  and  unity  gain.  In  other  words, 

the  partial  sum  f (t)  can  be  considered  as  the  output  of  the  Ideal  low 
n 

pass  filter  excited  by  f(t).  This  is  a more  explicit  Interpretation  of 
the  partial  sum  of  the  Fourier  series  and  Is  Important  for  deriving 
the  Gibbs'  phenomenon  and  windowing  effects  later  on. 


2.2  GIBBS'  PHENOMENON  IN  THE  TIME  DOMAIN 

Since  Gibbs  and  BScher^^,  many  Investigators  have  worked  on  the 
property  of  the  Fourier  series  of  a function  having  points  of  discon- 
tinuity. Here,  we  adopt  the  derivation  given  by  Gui 1 lemin^'^ . 

If  the  given  function  f(t)  Is  discontinuous  at  t ■ a,  the  property 
of  the  truncated  Fourier  series  In  the  vicinity  of  t * a is  equivalent 
to  that  of  replacing  f(t)  by  the  biased  step  function; 

f#(t)  ■ Au(t-a)  + f(a")  (17) 

where 

f (a  ) ■ 1 Im  f (a-e)  (18) 

e-*0 

and  A is  a real  constant,  such  that 

A * llm  [f(a+e)  - f(a-e)]  (19) 

e-*0 
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Considering  Eq.  (17);  the  second  term  f(a  ) Is  a constant  and  continuous; 
therefore,  we  can  drop  it  for  the  time  being  when  studying  the  proper- 
ties of  the  partial  sum  of  the  Fourier  series  for  a function  at  the 
point  of  discontinuity.  To  drop  the  term  f(a  ),  let 

g(t)  - f(t)  - f(a')  (20) 

For  n very  large,  It  Is  easily  shown  that 

g (t)  - fJO  ' *(•">  (*«) 

n n 


Now,  from  Eq.  (14),  we  have 


where 


a (t)  - f w (t-x)  g (t)  dt 
n n a 


g (t)  - f (t)  - f(a  ) - Au(t-a) 

o a 


<*>•/ 


(22) 


(23) 


Without  loss  of  generality,  let  us  assume  a ■ 0.  Thus  we  have 

Au(x)  dx  (24) 


According  to  the  properties  of  the  convolution  Integral,  we  have 

(x)  dx  (25) 


gn(t)  - A J wfl 


Substituting  Eq.  (15)  Into  Eq.  (25),  yields 


9 


sin  (nu.  x) 

-00 

f"V  , 

-i  iiiXdy 

* J-.  y 


In  terms  of  the  sine-integral  function,  Sl(x),  defined  as 


Si  (x) 


Eq.  (26)  can  be  rewritten  as 


9n(t)  “ A[l  + 7SI(nV}] 


The  above  function  is  shown  in  Figure  2-1.  Therefore,  as  n becomes 

very  large,  g (t)  has  a global  maximum,  occurring  at  t * ——  , 
n nco- 


9n (nwT^  * A[t+7S,(it)] 


- 1.0895  A 


and  • global  minimum,  occurring  at  t ■ - , 


-0.0895  A 


As  n-*»,  -2—  and  - -^—approach  t ■ 0 which  is  the  point  of  discontinuity 

nu»0  n«0 

and  the  ripples  In  Figure  2-1  are  compressed  into  a single  vertical  line 
at  t • 0.  Eqs.  (30)  and  (31)  correspond  to  the  top  and  bottom  of  this 
vertical  line  segment.  In  other  words,  the  overshoot  at  the  point  of 
discontinuity  Is  given  by 

1 im  | Max  g (t)  -A|  - 0.0895  A - 3%  A (32) 

n 

n-*» 

and  the  undershoot  Is  given  by 

11m  |Mln  gn(t)|  - 0.0095  A ~ 9%  A (33) 

n-*» 

The  Fourier  series  of  a discontinuous  function  results  In  a 9$  overshoot 
In  the  vicinity  of  the  point  of  discontinuity.  This  Is  known  as  the 
Gibbs'  phenomenon. 
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This  behavior  Is  Important  In  dealing  with  time  domain  analysis  of 
periodic  signals.  However,  In  the  study  of  digital  filters  or  discrete 
systems  the  same  phenomenon  Is  encountered  In  the  frequency  domain;  and 
It  Is  helpful  to  analyze  digital  filters  from  Gibbs'  viewpoint. 

2.3  GIBBS'  PHENOMENON  IN  THE  FREQUENCY  DOMAIN  - REAL  PERIODIC  FUNCTIONS 

Let  F (ui)  be  a real  periodic  function  of  frequency  u>  with  period  P 
P P 

and  bounded  In  (-  y,  y) . As  with  the  periodic  function  in  the  time 
domain,  we  can  expand  F(w)  as  a Fourier  series  In  the  frequency  domain: 


where 


and 


F(w) 


e‘Jkt0“ 


(34) 


JktQx 


dx 


(35) 


t 


0 


2ir 


(36) 


The  nth  partial  sum  is  defined  as 


FnU)  - 


”Jktow 

0 


(37) 


Let  us  use  the  second  Interpretation  of  the  partial  sum  of  the  Fourier 
series.  If  f(t)  Is  thi>  Inverse  Fourier  Transform  of  F(u>),  then 

f(t)  - ~ j F(w)  eJh)t  dm  (38) 
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F roii  Eq . ( 3*0  , we  have 


Mt)  .■  T.  <5  ( t-ktQ) 
k--“ 


(39) 


From  Eq.  (37).  the  inverse  Fourier  Transform  of  the  nth  partial  sum  is 


fn«> 


£ Ak  5<t-k,‘ 


m 


Now,  according  to  the  definition  of  the  Fourier  Transform, 


FnU) 


/•>. 


(t)  e'jwt  dt 


(41) 


r"tc 

J-nt, 


f(t)  e'jut  dt 


(42) 


and  substituting  Eq.  (38)  into  Eq.  (42),  using  the  dummy  variable  x for 
u in  Eq.  (38) , yields 


F (u) 

n 


1 r / ° f fM  •Jxt  ••J“t  dx  dt 

J-ntn  -*> 


. y f ,w  n .^>  dt  dx 
-»  "'•nt. 


(43) 


Eq.  (43)  can  be  rewritten  in  the  form  of  a complex  convolution 


F„<“> 


j_r 

2e  / 


F(x)  W (u-x)  dx 
n 


(44) 


13 


where 


2 sin  (ntotp) 
10 


(45) 


(46) 


may  be  called  the  aperiodic  Dirichlet  kernel  in  the  frequency  domain. 

Taking  the  inverse  Fourier  Transform  of  W (to),  yields 

n 


wn(t) 


(’  llf  nt0 

1 1 1 > ntn 


(47) 


According  to  the  complex  convolution  [Eq . (44)],  we  have 


f (t)  - f(t)  w (t) 
n n 

|f(t)  |t|<ntQ 

l0  1 1 1 > ntQ  (48) 

In  other  words,  wn(t)  ' s a rectangular  window  in  the  frequency  domain. 

Now,  assume  that  F(u)  Is  a discontinuous  function  with  a point  of 
discontinuity  at  id  * a.  Following  the  derivation  In  Section  2.2, 

F (u)  converges  to  the  biased  step  function  In  the  vicinity  of  id  ■ a, 

Fa(u>)  - Au(to-a)  + F(a")  (49) 

where 

F(a  ) • I Im  F(a-e) 
c-*0 
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(50) 


and  A ts  the  real  constant. 


A * lim[F(a+e)  - F(a-t)]  (50 

e-*0 

Since  convergence  In  the  vicinity  of  w ■ a is  governed  by  the  first  term 
of  Eg.  (49) , we  can  assume,  for  simplicity  and  without  loss  of  generality, 
that  a * 0 and  F(a”)  - 0.  Under  this  assumption  we  have 

F (w)  * Au(w)  (52) 

a 

From  Eg.  (44),  the  partial  sum  Is  given  by 


Fn(u») 


F (x)  W (urx)  dx 
a n 


u(x)  W (w-x)  dx 

n 


Substituting  Eg.  (46)  into  Eg.  (53),  yields 


(53) 


F («) 

n 


. f sin  [n 

. L u(x)  

U I U) 


tnt0(u-x)] 


dx 


. (u  tin  (nt0  y) 

•?  I — r2— 

••• 


or 


F (w) 


i rvup 

* •• 


dz 


(54) 


(55) 
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Comparing  Eq.  (55)  with  Eq.  (26),  we  find  that  convergence  of  a real 
periodic  function  in  the  frequency  domain  in  the  vicinity  of  the  point 
of  discontinuity  is  similar  to  that  of  a real  periodic  function  in 
the  time  domain.  Eq.  (55)  can  be  written  as 

Fn(w)  “ A[l+  7S,(nt0u)]  (56) 

Following  the  reasoning  in  Section  2.2,  the  Fourier  series  of  a real 
discontinuous  function  in  frequency  domain  results  in  a 9%  overshoot 
in  the  vicinity  of  the  point  of  discontinuity.  Therefore,  we  have 
determined  the  counterpart  of  the  Gibbs'  phenomenon  in  the  time  domain; 
that  Is,  the  Gibbs'  phenomenon  in  the  frequency  domain. 

We  now  have  identical  properties  for  the  truncated  Fourier  series 
of  real  periodic  functions  In  the  time  domain  and  the  frequency  domain: 

(1)  The  partial  sum  of  a real  periodic  function  in  the  time 
domain  can  be  expressed  as  the  convolution  of  the  original 
function  and  an  Ideal  filter  with  unity  gain;  the  partial 
sum  of  a real  periodic  function  in  frequency  domain  can  be 
expressed  as  the  complex  convolution  of  the  original  func- 
tion and  a rectangular  window  with  unity  height. 

(2)  The  Fourier  series  for  both  the  real  periodic  function  In 
the  time  domain  and  the  real  periodic  function  In  the  fre- 
quency domain  has  a 9%  overshoot  in  the  vicinity  of  the 
points  of  discontinuity;  In  other  words,  the  Gibbs' 
phenomenon  for  a real  periodic  function  in  the  frequency 
domain  is  exactly  the  same  as  that  for  a real  periodic  func- 
tion in  the  time  domain. 
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It  Is  Important  to  point  out  that  we  have  limited  our  study  to  the 
partial  sum  of  the  truncated  Fourier  series  for  real  periodic  functions 
In  the  frequency  domain.  In  reality,  we  often  must  consider  complex 
periodic  functions  in  the  frequency  domain.  For  this  case,  the  first 
result  listed  above  is  valid;  however,  Gibbs'  phenomenon  for  a complex 
periodic  function  in  the  frequency  domain  must  be  redefined. 

2.1*  GIBBS'  PHENOMENON  IN  THE  FREQUENCY  DOMAIN  - COMPLEX  PERIODIC 
FUNCTIONS 


Since  we  did  not  place  any  particular  restrictions  on  F(u>)  when  we 
derived  Eq.  (44)  for  the  partial  sum  of  a periodic  function  in  the  fre- 
quency domain,  Eqs.  (44)  and  (46)  are  still  valid  when  F(w)  is  a com- 
plex periodic  function.  However,  if  the  Gibbs'  phenomenon  is  to  be 
considered,  we  need  to  treat  the  real  part  and  imaginary  parts  indivi- 
dually. Let 


F (o»)  ■ FR(u)  + Jf  | (w) 


(57) 


The  window  function  in  Eq.  (46)  is  real;  therefore,  the  partial  sums  of 
Fr(u)  and  F( (u)  are 


FRn(w) 


"k  L 


Fr(x)  Wn(w-x)  dx 


(58) 


and 


n(w)  " k j 
«/-00 


F.(x)  W (w-x)  dx 
n 


(59) 


respectively.  Now,  If  F(u)  has  a point  of  discontinuity  at  u>  - a and 
letting 


C - ilm  (F(a+c)  - F(a-e)] 
e-*0 


(60) 
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wh«re  C Is,  In  general,  complex  and  can  be  defined  as 

C-CR+JC,  (61) 

where 


CR  - llm  [FR(a+e)  - F (a-e) ] 
e-*0 


and 


C,  - llm  (F  (a+c)  - F (a-c)] 
e-*0 


W 


(62) 


(63) 


Now,  we  can  treat  FR(w)  and  F((ui)  as  two  real  periodic  fun«lons  In 
the  frequency  domain.  From  Eq.  (55)  we  have  the  partial  sums 


c.  r'o"  , 

F.J“>  * T-  I «*> 

" “00 

and 

c,  rnto“  , 

Fln(w)"7-|  1Tidr  (65) 

" -00 

In  other  words.  In  the  vicinity  of  a point  of  discontinuity,  the  partial 
sum  of  the  complex  function  In  the  frequency  domain  is 


Fn(«)  - FRf)  (w)  + J F | n (*°) 


18 


Eq.  (66)  can  be  rewritten  as 

FnU)  ’ C[l+  7*,("ty*)]  (67) 

Therefore,  the  Fourier  series  for  the  complex  periodic  function  In  the 
frequency  domain  results  in  an  overshoot,  CR  £|-SI(Tr)-lJ  - 9%  CR,  in 

the  real  part  and  an  overshoot,  C|£^-  Si  (ir)-lj  ■ 3%  Cj,  in  the  Imaginary 
part  near  the  point  of  discontinuity.  Although  the  magnitudes  of  the 
overshoots  in  the  real  and  Imaginary  parts  are  generally  not  equal,  the 
percentages  are  still  the  same  as  the  Gibbs'  phenomenon  In  the  time 
domain  if  CR  » C( . 

When  dealing  with  the  transfer  function  of  digital  filters  or  dis- 
crete-time control  systems,  we  are  often  interested  in  the  magnitude 
and  phase  relationships  of  the  frequency  response  instead  of  the  real 
part  and  the  Imaginary  part.  In  terms  of  phase  and  magnitude,  the 
Gibbs'  phenomenon  Is  more  complicated.  Figure  2-2  shows  the  situation 
for  overshoots  of  magnitude  and  phase  relative  to  those  of  the  real  part 
and  the  Imaginary  part.  The  Jump-amplitude  at  the  point  of  discontin- 
uity is  given  by 


AM  - | F (a+) | - | F (•”) | 

The  jump-phase  Is 


A9  - ARG[F(a+)l  - ARG[F(a')] 
The  overshoot  in  amplitude  Is 

«,  - |F(a+)  + 9*C|  - |F(a+)l 
and  the  undershoot  In  amplitude  is 

«2  - |F(a‘)f  - |F(a")  - 9«| 


(68) 


(69) 


(70) 


(7D 
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The  overshoot  in  phase  Is 


0,  - ARG[F (a+)  + 9«]  - ARG[F(a+)  ] (72) 

and  the  undershoot  in  phase  is 

02  - ARG[F(a") ] - ARG[F (a")  - 9%C]  (73) 


FIGURE  2-2.  OVERSHOOTS  OF  MAGNITUDE  AND  PHASE  RELATIVE  TO 
THOSE  OF  REAL  AND  IMAGINARY  PARTS 

61 

In  general  6j  i fij,  0^  4 the  percentages  of  overshoot  (^j)  and  under- 
^2 

shoot  (^yj)  in  amplitude  are  not  3%  except  in  some  special  cases. 

2.5  DEMONSTRATION  OF  GIBBS'  PHENOMENA  WITH  HP  FOURIER  ANALYZER 

In  this  section,  Gibbs'  phenomena  in  the  time  domain  and  the  fre- 
quency domain  will  be  demonstrated  through  the  use  of  the  Hewlett  Packard 
Fourier  Analyzer.  The  Fourier  Analyzer  performs  the  FFT  algori thm^  to 
compute  the  Discrete  Fourier  Transform  (DFT)  of  discrete  signals  in  the 
time  domain  or  in  the  frequency  domain.  The  DFT  can  be  used  to  evaluate 
the  coefficients  of  Fourier  series  from  the  periodic  waveform  or  to 
synthesize  the  waveform  from  the  coefficients  of  the  Fourier  series 
(Appendix  A) . 
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To  demonstrate  the  Gibbs'  phenomenon  In  the  time  domain,  the  following 
wave  was  chosen  as  an  example: 

!1  -y  < t < 0 

0 1 “ °*  T (7A) 

-1  0 < t < i 


The  square  wave  defined  by  Eq.  (7*0  is  characterized  by  an  amplitude  of 
unity  and  a period  of  I second.  Program  No.  I (Appendix  C)  Is  a listing 
of  the  program  utilized  for  this  demonstration.  The  number  of  samples 
required  for  the  FFT  algorithm  was  **096.  Figure  2-3  shows  the  resulting 
waveforms  for  various  partial  sums  of  the  Fourier  series  for  Eq.  (7**). 
The  maximum  overshoots  for  different  nth  partial  sums  are  shown  in  Fig- 
ure 2-*i.  The  maximum  overshoot  approaches  0.181  when  n Is  very  large. 
The  height  of  the  jump  at  the  point  of  discontinuity  Is  2;  thus,  the 
ratio  of  overshoot  to  the  height  of  the  jump  Is  0.0905  or  9.05%.  This 
is  almost  equal  to  the  prediction  of  8.95%  based  on  Gibbs'  phenomenon. 


The  following  example  Is  presented  to  demonstrate  the  Gibbs' 
phenomenon  In  the  frequency  domain.  Consider  the  periodic  function 
F(«)  » FRL)  ♦ JF(((i>),  where 


Fr(w)  " 


1.0  M < 0.25 


-1.0  0.25  < M < 0.5 


/-2.0  + 4.0a)  -0.5  < u>  < -0.25 


F j (to)  ■ \ **.0w 


LI  < 0.25 


2.0  - *i.0(i)  0.25  < w < 0.5 
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The  period  Is  I radian.  The  Fourier  Analyzer  only  displays  FR(w)  and 
Fj  (w)  for  0 <^0.5  and  assumes  FR(-w)  ■ FR(w),  Fj  (-w)  * Fj(w).  Fr(m) 

and  F|(u)  as  well  as  the  magnitude  and  phase  of  F(w)  are  shown  In 
Figures  2-5 (a)  and  (b) . The  results  of  various  partial  sums  are  shown 
in  Figures  2-5(c)  through  2-5 (g) - The  magnitude  of  F(u>)  displays  jumping- 
up  and  jumping-down  spikes  at  the  middle  of  the  curve.  The  jumping-down 

spike  is  due  to  the  real  part  of  the  Fourier  series  [Fd(u)]  converging 

1 * I 

to  0 at  t * jj-,  and  the  magnitude  of  F(u)  converging  to  1.0  at  t • jp. 

The  jump-up  spike  is  due  to  Gibbs'  phenomenon.  The  maxi mums  of  the 

ripple  magnitude  for  the  jump-up  spike  of  different  nth  partial  sums 

are  shown  in  Figure  2-6;  the  peaks  of  the  ripples  approach  1.542.  In 

other  words,  the  ratio  of  the  maximum  of  the  ripples  to  the  jumping  height 

in  the  real  part  FR(w)  (which  Is  2.00)  Is  1(1. 5422  - 1.0002),/2  - I]  * 

2.00  ■ 0.087  or  8.7$.  This  is  very  close  to  the  prediction  of  Gibbs' 

phenomenon.  Since  F( (w)  is  continuous,  the  overshoots  at  t • are  due 

to  the  discontinuity  of  FR(u>)  only.  The  magnitude  function,  as  shown 

in  Figure  2~5(b),  Is  also  continuous;  therefore,  taking  the  ratio  of 

overshoot  to  the  jumping  height  (which  is  0)  of  the  magnitude  function, 

we  find  that  the  ratio  is  infinite.  The  result  indicates  that  the  ratio 

of  overshoot  to  the  jumping  height  in  the  magnitude  function  for  a 

complex  periodic  function  in  the  frequency  domain  is  not  always  equal 

to  3%. 

2.6  REHARKS 

Starting  with  the  truncated  Fourier  series,  we  have  defined,  derived, 
and  interpreted  Gibbs'  phenomenon  in  the  time  domain  and  frequency 
domain.  The  latter  Is  a counterpart  of  the  former.  A one-to-one 
correspondence  exists  with  the  exception  that  Gibbs'  phenomenon  in  the 
frequency  domain  deals  with  complex  periodic  functions  while  in  the 
time  domain  Gibbs'  phenomenon  deals  with  real  functions. 

This  chapter  offers  a basis  for  nonrecursive  filter  analysis  via  the 
windowing  method. 
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(f)  128th  PARTIAL  SUM 

FIGURE  2-5.  THE  PERIODIC  FUNCTION  F(u)  AND  THE  MAGNITUDE  ANO  PHASE  OF 
VARIOUS  PARTIAL  SUMS  OF  F(w)  (Continued) 
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FIGURE  2-6.  THE  JUMP-UP  SPIKE  OF  DIFFERENT  PARTIAL  SUMS  OF  EXAMPLE  2 
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CHAPTER  III 

THE  GENERALIZED  GIBBS'  PHENOMENON  AND  WINDOWS 

3.0  INTRODUCTION 

In  Chapter  1 1 , we  have  shown  that  the  truncation  of  a Fourier  series 
can  be  considered  the  same  as  the  windowing  of  a function  through  the 
use  of  a rectangular  window.  If  a point  of  discontinuity  occurs  in  the 
function,  then  Gibbs'  phenomenon  results.  To  eliminate  or  reduce  the 
Gibbs'  phenomena,  rectangular  windows  must  be  modified  into  windows  of 
other  shapes.  When  modified  windows  are  used,  the  corresponding  effects 
are  referred  to  as  the  generalized  Gibbs'  phenomena.  In  turn,  based 
on  the  generalized  Gibbs'  phenomena,  we  will  develop,  classify,  and 
Investigate  various  windows  both  in  the  time  domain  and  frequency  domain; 
in  particular,  we  will  establish  a trapezoidal  window.  It  is  believed 
that  this  Is  the  first  time  that  the  trapezoidal  window  has  been  uti- 
lized to  offer  a unified  picture  of  windowing  techniques. 

3.1  WINDOWING  IN  THE  TIME  DOMAIN 

Let  us  restate  the  nth  partial  sum  of  the  truncated  Fourier  series 


f (t)  - I f(x)  w (t-x)  dx 
n I _ n 


- f(t)  * w (t) 
n 


where  Is  real  convolution  and 


wn(t)  * 


sin  (nw.t) 


1 


Therefore,  since  It  Is  apparent  that  Eq.  (2)  represents  an  ideal  filter 

with  a cutoff  frequency  nu>g,  we  can  say  that  the  partial  sum  f n ( t)  is 

governed  by  the  characteristics  of  the  ideal  filter  w (t).  As  shown  in 

n 

figure  3-1  (a),  wn(t)  displays  low-damped  oscillations  around  t ■ 0;  result 

ing  in  (1)  a slow  convergence  rate  and  (2)  Gibbs'  phenomenon.  To  improve 

the  convergence  rate  and  eliminate  the  Gibbs'  phenomenon,  we  must  use 

other  filters  or  window  functions  Instead  of  wn(t)  given  In  Eq.  (2). 

If  we  filter  the  nth  partial  sum  f (t)  with  a filter  with  Impulse 

n 

response  wtn(t)»  the  resulting  smoothed  waveform  is  given  by 

f*(t)  - f (t)  * w (t)  (3) 

n n tn 

Vv. 

Substituting  Eq.  (l)  into  Eq.  (3),  yields 

f*(t)  " f(t)  * w (t)  * w (t)  (U) 

n n tn 

Since  the  convolution  integral  is  an  associative  operation,  Eq.  (it)  can 
be  rewritten  as 

f*(t)  - f(t)  * w(t)  (5) 

where 

w(t)  - w (t)  * w (t)  (6) 

n tn 

Let  us  assume  that  the  Fourier  Transform  of  w.(t)  is  W (w)  and 

2 tn  tn 

Wtn(u)  ■ 0 for  u > nwg,  where  Uq  • T •*  the  period  of  f(t). 

Taking  the  Fourier  Transform  of  w(t) , we  have 

W(w)  ^ Flw(t)J  • V (w)  W («)  (7) 

n tn 
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Since  (u»)  represents  an  ideal  filter  with  cutoff  frequency  nuig, 

W(u>)  - Wtn(u>)  (8) 

and, 

w(t)  " wtn  ^ (9) 

Eq.  (5)  now  becomes 

f*(t)  * f(t)  * wtn(t)  (10) 

The  result,  f^n (t) , may  be  called  the  windowed  partial  sum  of  the  periodic 
function  f(t).  From  this  viewpoint,  truncation  of  the  Fourier  series 
(Chapter  li)  can  be  considered  as  a special  case  of  windowing. 

3.2  THE  GENERALIZED  GIBBS  PHENOMENON  IN  THE  TIME  DOMAIN 

As  pointed  out  earlier,  the  property  of  the  truncated  Fourier 
series  in  the  vicinity  of  a discontinuous  point  is  equivalent  to  that 
of  replacing  f(t)  by  a biased  step  function.  If  only  the  oscillations 
in  the  vicinity  of  the  point  of  discontinuity  are  considered,  we  can 
replace  f(t)  by 

f#(t)  - Au(t)  (II) 

where  A Is  the  jump  height  of  f(t)  at  the  discontinuous  point  t > 0. 

When  the  periodic  function  f(t)  Is  windowed  by  wtr|(t)»  we  can  apply 
this  Idea  and  evaluate  the  overshoot  in  the  vicinity  of  the  discon- 
tinuous point.  Let  us  define 

9n(t)  - f#(t)  * wtn (t)  (12) 
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r 


Then 


gn(t)  - A J u(x)  wtn 

w •«* 


(t-x)  dx 


(13) 


•*£ 


A I w (x)  dx 
tn 


om 


To  facllltete  our  discussion,  let  us  define  the  unit  window  w^(t)  for 

the  class  of  windows  w (t) , for  n » I,  2,  ....  such  that 

tn 


Vt(w)  * F[wt(t)]  - 0 for  u » 1 


(15) 


and 


Thus,  we  have 


V“>  2 ’ W«(^) 


w (t)  - no>0  wt(nw0t) 


Substituting  Eq.  (17)  Into  Eq.  (1*0  yields 


Letting 


3 (t)  - A f nu»0  wt(nu0x)  dx 
J —00 

fnV 

wt(y)  dy 

-00 

g(t)  - A j w{(y)  dy 
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(16) 


(17) 


(18) 


(19) 


9n(t)  - g(nw0t)  (20) 

Assume  that  g(t)  has  several  ripples  with  the  local  maxima  occurring  at 
t * tj,  t^,  ....  then  gn(t)  has  local  maxima  occurring  at 

t ■ — — , — — , ....  Therefore,  when  n-*«,  all  the  ripples  of  g (t)  are 
nw0  0 n 
compressed  into  a single  vertical  line  at  t - 0.  This  may  cause  the 

same  effect  as  Gibbs'  phenomenon. 

t 

In  order  to  obtain  the  same  jump  height  (A)  after  windowing  for 
n-M°,  we  need  to  have 


11m  [llm  g (e) ] • A f w (y)  dy  - 
e-K)+  n—  n 1 


Wt(0)  - J wt(y)  dy  - 1 


In  other  words,  the  necessary  conditions  for  an  acceptable  time  domain 
window  are  the  following: 


(1)  Its  Fourier  Transform  Wtn(o>)  vanishes  for  |u|  > nuQ. 


(2)  W.  («)  - I for  w - 0. 
in 


(3)  • “,„(■) 


Let  us  define 


g*(t) 


i 


and 


9C  * l.u.b.  (g*(t)}  - 1 t € (0,-)  (2*») 

Then 

gc  - J l.u.b.  (gn(t)}  - I t CIO,-)  (25) 

From  Eqs.  (19),  (23)  and  (24),  and  fol lowing  the  same  arguments  used  In 
Section  2.1  of  Chapter  1 1 , we  can  define  the  generalized  Gibbs'  phenom- 
enon In  the  time  domain  as: 

(1)  If  a periodic  function  has  a discontinuous  point  with  jump 
height  (A)  and  Is  windowed  by  a time  domain  window  Mtn(t), 
then  the  overshoot  in  the  vicinity  of  the  point  of  discon- 
tinuity is  given  by 

t » 9 A (26) 

c 

(2)  The  ratio  of  the  overshoot  to  the  Jump  height  is  a con- 
stant, 3C,  which  is  dependent  on  the  unit  window  wt(t) 
for  the  specific  class  of  windows  generated  by  wtn(t)  - 
wt(nwQt). 

The  unit  window  for  the  rectangular  window  or  the  ideal  filter  is 

( 1 l“l  i 1 

Wt(«)  - j (27) 

( 0 |u|  > 1 

and 

wt(.)  - ii=-L  (26) 
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1 .u.b.  {g*(t) } ■ 1 im  g*(t) 

t-*» 


j wf(t)  dt 


In  other  words  If  wt(t)  or  wtn(t)  is  non-negative,  g£  ■ 0. 

Consider  the  triangular  window  which  has  the  unit  window 


(33) 


Wt(w)  - 


|-u  0 <_  (i)  <_  1 

l+w  -I  < u>  < 0 


(3A) 


The  corresponding  time  function  is  given  by 


*t(t)  ■ 


sin 


li) 


t 

2 


(35) 


It  Is  evident  that  w^(t)  >_ 0 for  t G (0,»).  Therefore  g£  ■ 0 and  we 
may  conclude  that  the  triangular  window  can  be  used  to  eliminate  the 
Gibbs'  phenomenon  completely. 


3.3  THE  TRAPEZOIDAL  WINDOW  IN  THE  TIME  DOMAIN 


Consider  the  class  of  trapezoidal  windows  with  unit  window 


wtU) 


where  0 < b < 0.5 


I - (i) 

ir~ 

i 

I ♦ to) 

“5E 


1 - 2b  < u <_  1 
|u|  <.  1 - 2b 
-1  < id  < -(1  - 2b) 


(36) 
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l 


(37) 


The  corresponding  time  function  is  given  by 


M 

fs In  (bt)l 

fs In  ((1  - b)  t)l 

L « J 

L bt  J 

L (1  - b) t J 

From  Eqs.  (36)  and  (37)  we  recognize  that  the  trapezoidal  windows 
become  rectangular  windows  If  b - 0 and  triangular  windows  if  b « 0.5. 
Therefore,  the  trapezoidal  window  Is  a general Izatlon  for  straight  seg- 
ment windows.  The  overshoot  ratio  is  gc  - 8.95?  when  b - 0 and  g - 0 
when  b - 0.5.  Figures  3-1  through  3-6  show  the  trapezoidal  window 

wt(t)  and  J f w£(c)  dc  for  various  b's.  Program  No.  3 (Appendix  C)  was 

used  for  this  purpose.  The  overshoot  ratios  g£  for  different  b's  are 
shown  In  Table  3“ 1 and  Figure  3”7. 


The  width  of  the  main  lobe  is  defined  as  twice  the  Interval  from 
t - 0 to  the  first  zero-crossing  of  wt(t).  From  Eq.  (37)  we  find  that 
the  width  of  the  main  lobe  is  given  by 


(38) 


for  the  unit  window.  Figure  3-8  shows  T*  as  the  function  of  b.  If  the 
trapezoidal  window  has  a cutoff  frequency  n«0,  the  width  of  the  main  lobe 
is  given  by 


(39) 


From  Figure  3-1  through  3-6  we  see  that  as  b Increases,  not  only 
does  the  overshoot  ratio  gc  become  smaller  and  smaller,  but  also  the 
side  lobe  becomes  smaller  and  damps  out  more  quickly.  From  this  view- 
point, the  triangular  window  is  the  best  window  for  the  overshoot  ratio 
gc  - 0 since  the  side  lobe  damps  out  quickly.  However,  If  we  examine 
Figure  3-8,  we  find  that  the  width  of  main  lobe  Increases  with  b;  there- 
fore, the  width  of  the  main  lobe  Is  largest  for  the  triangular  window. 
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FIGURE  3*1*.  TIME  AND  STEP  RESPONSES  OF  TRAPEZOIDAL  WINDOW  WITH  b 


FIGURE  3-6.  TIME  AMD  STEP  RESPONSES  OF  TRAPEZOIDAL  WINDOW  WITH  b 


FIGURE  3-7.  OVERSHOOT  RATIO  OF  THE  TRAPEZOIDAL  WINDOW  WITH  RESPECT  TO 

THE  PARAMETER  b 


FIGURE  3-8.  WIDTH  OF  THE  HA IN  LOBE  OF  THE  TRAPEZOIDAL  WINDOW  WITH 

RESPECT  TO  THE  PARAMETER  b 


r 


From  the  properties  of  the  convolution  Integral,  an  increase  in  the 
width  of  the  main  lobe  means  a decrease  In  the  resolution  of  the  windowed 
waveform.  A good  window  is  characterized  by  the  following:  (1)  the 

width  of  the  main  lobe  Is  small,  (2)  the  side  lobes  damp  out  quickly 
and  (3)  the  overshoot  ratio  is  small  or  zero. 

TABLE  3-1 

OVERSHOOT  RATIO  g£  AND  THE  WIOTH  OF 
MAIN  LOBE  T* 


b 

9C(S) 

T* 

w 

271 

0 

9.0 

1.00 

0.05 

8.8 

1.05 

0.10 

8.7 

1.11 

0.15 

8.4 

1.18 

0.20 

7.9 

1.25 

0.25 

7.0 

1.33 

0.30 

5.9 

1.43 

0.35 

4.3 

1.54 

0.40 

2.1 

1.67 

0.45 

0.7 

1.82 

0.50 

0.0 

2.00 

Consider  the  square  wave  discussed  in  Chapter  II.  The  partial 
sum  converges  slowly  and  the  Gibbs'  phenomenon  with  an  overshoot  ratio 
of  approximately  9%  occurs.  Figure  3”9(a)  shows  the  15th  partial  sum 
of  Its  truncated  Fourier  series.  Figures  3-9(b),  (c)  and  (d)  show  the 
15th  partial  sums  of  the  windowed  Fourier  series  utilizing  trapezoidal 
windows  with  b ■ 0.33*  0.4  and  0.5  respectively.  From  Figure  3"9  we 
see  that  the  underdamped  oscillation  Is  eliminated  and  the  overshoot 
spike  is  reduced  after  windowing  by  trapezoidal  window  with  b - 0.33. 
Figures  3-10(a)  through  3“ 10 (d)  show  the  255th  partial  sums  of  the 

I 

£ 

» 

At 

I 


windowed  Fourier  series  utilizing  trapezoidal  windows  with  b ■ 0,  0.33, 
0.4  and  0.5  respectively.  From  Figures  3~3  and  3-10,  the  effect  of 
windowing  is  evident. 

3.4  WINDOWING  AND  THE  GENERALIZED  GIBBS'  PHENOMENON  IN  THE  FREQUENCY 
DOMAIN 

In  the  frequency  domain  we  have  a complex  periodic  function  to 
deal  with  and  we  can  develop  windowing  techniques  accordingly.  Con- 
sider Eq.  (44)  In  Chapter  II,  rewritten  as 

F (id)  • F(id)  ® W (id) 
n n 

where  " ® " implies  complex  convolution  and 

2 sin  (nwtA) 

W>) = — <*0 

n b) 

is  the  Dlrichlet  kernel  or  rectangular  window  in  the  frequency  domain. 
Following  the  same  reasoning  given  In  Section  3-1,  If  the  nth  partial 
sum  F (w)  is  windowed  by  a frequency  domain  window  V^(u)  with 
w^  (t)  having  the  nonzero  interval  |t|  <_ ntQ,  it  yields 

F*(io)  • pn(“)  ® 

■ F (cd)  ® W^n(fa>)  (42) 

In  order  to  study  the  oscillation  In  the  vicinity  of  the  point  of 
discontinuity  after  windowing  by  Wfn(w),  we  replace  F(w)  by 

F (u)  ■ Cu  (to)  (43) 

c 
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FIGURE  3-9.  EFFECT  OF  TRAPEZOIDAL  WINDOWING  ON  THE  15th  PARTIAL  SUM  OF  THE  TRUNCATED  FOURIER  SERIES 


EFFECT  OF  TRAPEZOIDAL  WINDOWING  ON  THE  255th  PARTIAL  SUM  OF  THE  TRUNCATED  FOURIER  SERIES 


where  C - CR  ♦ Cjj,  end  CR,  C(  are  the  jump  heights  of  FR(w)  and  F(  (u>) 
at  the  discontinuous  point,  respectively;  and  u(u)  is  the  unit  step 
function  of  u.  Let  us  define 


GnU)  - FcU)  ® Wfn(u>) 


u (x)  Wfn(u-x)  dx 
Wfn(x)  dx 


(kk) 


Similar  to  the  time  domain  window,  we  can  define  the  unit  window  W^(w) 
for  the  class  of  windows  W^n(u>),  n * I,  2,  ...,  such  that 


w^(t)  - F~'[Wf(w)l  ■ 0 for  t > 1 (**5) 


and 


From  Eq.  (1*6)  we  have 


Wfn(u)  - ntQ  Wf(nt0u) 


Substituting  Eq.  (1*7)  into  Eq.  (M*)  we  have 


57  £ «o  Mf 


(ntQx)  dx 


c f'°“ 

■57  J 

-oa 
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(46) 


(47) 


Wf(y)  dy 


(48) 


If  we  let 


6n(w)  " h /:  Wf(y)  dy  (li 9) 

Then  we  have 


Gn(u)  - G(ntQai)  (50) 


If  G(«d)  has  ripples  with  local  maxima  occurring  at  u>  ■ o>j , 

id.  iD- 

then  G (id)  has  local  maxima  occurring  at  id  - -7—  , , ....  Therefore, 

n ntQ  ntQ 

as  n-w,  all  the  ripples  of  Gn(w)  are  compressed  Into  a single  vertical 
line  at  id  * 0.  In  order  to  obtain  the  same  jump  height  (C)  after  win- 
dowing for  n-*»,  we  need  to  have 

Urn  11m  G (c)  ■ §—  / Wf(y)  dy 

e-0  n-*»  n ZlT  J-o*  f 

. - C (51) 


wf(o)  mh  f_  Vy)  dy  * 1 (52> 

Therefore,  the  necessary  conditions  for  an  acceptable  window  in  the 
frequency  domain  are  the  following: 

(1)  Its  Inverse  Fourier  Transform  w^n(t)  vanishes  for  |t|  > ntQ 

(2)  Wf„(t)  - 1 for  t - 0. 

(3)  wfn(-t)  - -wfn(t) 


1*6 


Let  us  defli 


6*(u)  • j f Wf(y)  dy  - 
J -ao 


G ■ l.u.b.  {G*(u)>  - 1 u C (0,«) 
c 


Then,  from  Eq.  (50),  we  have 


G • i-  l.u.b.  (G  (w)}  - 1 w C (0,«) 
c t n 


Me  can  new  define  the  generalized  Gibbs'  phenomenon  in  the  frequency 
domain; 


(1)  If  a complex  periodic  function  In  the  frequency  domain  has 
a point  of  discontinuity  with  Jump  heights  CR  In  the  real 
part  and  C(  In  the  Imaginary  part,  and  is  windowed  by  a 
frequency  domain  window  W^n(<u),  the  overshoot  In  the 
vicinity  of  the  point  of  discontinuity  Is  given  by 


AR  " 6cCR 


In  the  real  part. 


A.  ■ G C.  in  the  Imaginary  part. 

i C I 


(2)  The  ratio  of  overshoot  to  the  Jump  height,  defined  as 
A.  A, 

r and  7—,  Is  a constant  G which  Is  real  and  dependent  on 
n Cl 

the  unit  window  W.U)  for  the  class  of  windows  generated  by 


Wfn (w)  • Wf(ntQu). 


The  unit  window  for  the  rectangular  window  is 


wf(t)  - 


Itl  « 1 


Itl  > 1 


Wf  («) 


2 sin  <u 


From  Eq.  (53)  we  have 


•*«■>  ■ h f“  * 


Gc  - l.u.b.  {G*(w) } - 1 


- TiilLXdy 


for  w C (0,«>) 


“ w Sl  (ir)  "I 


Gfi  in  Eq.  (61)  is  about  8.951.  This  Is  the  overshoot  ratio  for  the 
Gibbs'  phenomenon  In  the  frequency  domain  caused  by  truncating  the 
Fourier  series. 

In  order  to  eliminate  the  generalized  Gibbs'  phenomenon,  we  must 
choose  a window  such  that 


G ■ l.u.b.  (G*(w)}  - 1 » 0 
c 


i 


or 


l.u.b.  {G* (w) } - 1 


For  a class  of  windows  with  unit  window 


W^(m)  > 0 for  u C (0,«) 


than  G*(w)  Is  monotonlc  Increasing,  and  Gc  ■ 0.  Therefore,  Eq.  (63) 
gives  the  sufficient  condition  for  eliminating  the  generalized  Gibbs' 
phenomenon.  The  triangular  window  with  the  unit  window  defined  as 


wf(t)  - 


I - t 0 < t < 1 


1 + t -1  < t < 0 


has  a corresponding  frequency  function  given  by 


Wf(«)  - 


(?) 


Since  Wf(w)  >_  0 for  u C (0,»),  the  triangular  window  can  be  used  to 
eliminate  the  Gibbs'  phenomenon. 

When  overshoots  of  magnitude  and  phase  Instead  of  the  real  and 
Imaginary  parts  of  the  windowed  complex  period  functions  are  concerned, 
Eqs.  (70)  through  (73)  of  Chapter  II  are  valid  except  that  the  overshoot 
ratio  (91)  in  those  equations  must  be  replaced  by  G£  of  the  window  which 
Is  used. 

3.5  THE  TRAPEZOIDAL  WINDOW  IN  THE  FREQUENCY  00MAIN 


Consider  the  trapezoidal  window  In  the  frequency  domain  with  the 


*9 


unit  window  defined  as 


wf  (t) 


1 - t 
2b 

I 

1 + t 
“2b~ 


1 - 2b  < t < 1 
1 1 1 < 1 - 2b 
-1  < t < -(1  - 2b) 


(66) 


where  0 < b < 0.5.  The  Fourier  Transform  of  Eq.  (66)  Is  given  by 


wfM  ■ 2C-» 


(67) 


Figures  3~ 1 1 through  3“ 1 6 


show  the  amplitude  response  of 


Wf(a>) 

W^TST 


for 


w >,  0 and  for  b ■ 0,  0.1,  0.2,  0.3,  0.4  and  0.5,  respectively,  In  both 
the  linear  and  the  logarithmic  scale.  For  b - 0,  as  shown  In  Figure  3-11, 
the  side  lobes  are  underdamped;  as  a matter  of  fact,  It  Is  a rectangular 
window.  The  peak  of  the  highest  side  lobe,  denoted  by  S , is  -13.3  dB. 

As  b Increases,  S_  decreases.  When  b Is  0.418,  as  shown  In  Figure  3-17, 

m 

S_  approaches  a minimum  of  -32  dB.  When  b is  greater  than  0.418,  S 
*n  m 

increases.  When  b - 0.5,  a triangular  window  results  with  S - -26.5  dB. 

m 

Table  3~2  and  Figure  3-18  show  the  peak  of  the  highest  side  lobe  as  the 
function  of  parameter  b. 


The  bandwidth  of  the  main  lobe,  from  Eq.  (67),  Is  given  by 

Bw ' rH  (68) 

and  Is  the  same  as  T*  shown  In  Figure  3“8.  The  overshoot  ratio  G for 
the  generalized  Gibbs'  phenomenon  Is  the  same  as  g£  which  Is  the  over- 
shoot ratio  for  the  generalized  Gibbs'  phenomenon  In  the  time  domain  as 
shown  In  Figure  3~7.  From  Ftgures  3~7  and  3~18  we  see  that  for  0.38  < 
b < 0.5,  the  overshoot  ratio  G_  Is  less  than  33  and  S_  is  less  than 
-26  dB;  therefore,  the  trapezoidal  windows  with  0.38  < b < 0.5  give  good 
windowing  resulting  In  a small  overshoot  ratio  and  overdamped  side  lobes. 
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FIGURE  3-11.  TRAPEZOIDAL  WINDOW  IN  THE  FIGURE  3-12.  TRAPEZOIDAL  WINDOW  IN  THE 

FREQUENCY  DOMAIN,  b - 0 FREQUENCY  DOMAIN,  b = O.l 


TABLE  3-2 


PEAK  OF  THE  HIGHEST  SIDE  LOBE 

4 

b 

S (dB) 

fP 

0 

-13-3 

0.1 

-13.6 

0.2 

-15.1 

0.3 

-19.0 

0.35 

-22.7 

0.4 

-28.8 

0.418 

-32.0 

0.42 

-31.6 

0.45 

-28.4 

0.50 

-26.5 

From  Figure  3” 1 7 for  b - 0.418,  the  peak  of  the  highest  side  lobe 
Is  -32  dB  which  ts  the  lowest  Sffl  as  shown  In  Figure  3- 1 8 . This  Is  the 
optimal  trapezoidal  window  and  shall  be  used  In  the  design  of  nonrecur- 
slve  filters  and  compensators. 

Program  No.  A (Appendix  C)  was  used  to  determine  the  amplitude 
response  of  the  trapezoidal  windowed  Fourier  Transform  of  Eq.  (I). 

The  well  known  Hanning  and  Hamming  windows  are  defined  as: 

(I)  Unit  window  for  Hanning  window 

w(t)  ■ 0.5  + 0.5  cos  (irt)  1 1 1 < 1 

- 0 |t | > 1 (69) 


sin  m J_  sin  U+fa))  1 sin  U-tn) 

w 7 ir+w  7 ir-o) 


(2)  Unit  window  for  Hamming  window 


w(t)  - 0.54  + 0.46  cos  (rrt ) |t|  < I 

“ 0 |t  | > 1 (71) 

and 

W(w)  -1.08  + o.46  + 0.46  H-n.  .(V.m)-  (72) 

U 1T+(1)  1T-0) 

Figures  3-19  and  3“20  show  the  amplitude  response  of  the  Hanning  and 
the  Hamming  window  respectively.  The  peak  of  highest  side  lobes  of  the 
Hanning  window  Is  -31*5  dB  which  is  somewhat  greater  than  the  peak  of 
the  optimal  trapezoidal  window.  For  the  Hamming  window,  is  equal  to 
-43.2  dB  and  is  much  lower  than  the  Hanning  window.  However,  the  band- 
width of  the  main  lobe  is  3*44  it  for  the  optimal  trapezoidal  window, 
and  4ir  for  both  the  Hanning  and  the  Hamming  windows. 

Considering  the  second  example  at  the  end  of  Chapter  II  as  an 

illustration  of  the  effect  of  windowing  In  the  frequency  domain.  Figures 

3-21  and  3“25  show  the  magnitudes  and  phases  of  the  truncated  Fourier 

series  windowed  by  a rectangular  window  for  n - 16  and  256,  respectively. 

The  results  of  windowing  by  trapezoidal  windows  with  b » ~ and 

y for  n « 16  are  shown  in  Figures  3“22,  3-24  and  3"25  respectively. 

Figures  3~26,  3*27  and  3~28  show  the  results  of  windowing  by  trapezoidal 
171  212  1 

windows  with  b - jryy  , jsyy  and  y,  respectively,  for  n - 256.  Examina- 
tion of  these  waveforms  indicates  that  trapezoidal  windows  with  b > 0.3 
are  acceptable  for  windowing  in  the  frequency  domain. 

3.6  REMARKS 

Based  on  the  analysis  of  the  generalized  Gibbs'  phenomenon,  we  have 
explained  important  functions  of  various  windows  both  In  the  time  domain 
and  In  the  frequency  domain.  For  simplicity  In  exploring  the  idea,  the 
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trapezoidal  window  was  used  to  offer  a unified  picture.  In  the  following 
chapters,  we  shall  show  that  this  approach  is  powerful  and  helpful  in 
design  work. 


E 3-22.  THE  16th  PARTIAL  SUM  WITH  THE  TRAPEZOIDAL  WINDOW 
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FIGURE  3-24.  THE  16th  PARTIAL  SUM  WITH  THE  TRAPEZOIDAL  WINDOW 


FIGURE  3-28.  THE  256th  PARTIAL  SUM  WITH  THE  TRAPE20IDAL  WINDOW, 


CHAPTER  IV 

DESIGN  OF  NONRECURSIVE  FILTERS  USING  THE  WINDOW  METHOD 


4.0  INTRODUCTION 

Digital  filters  can  be  divided  Into  two  general  classes:  recursive 

and  nonrecursive  filters.  In  control  systems  applications,  the  non- 
recurslve  filter  has  two  important  advantages; 

1.  A nonrecursive  filter  has  a linear  phase  function  with 
poles  located  at  z * 0,  therefore,  a higher  degree  of 
stability  Is  inherent  In  the  system. 

2.  Since  the  Impulse  response  of  a nonrecurslver  filter  Is 
finite  In  length,  the  filter  has  only  a finite  memory,  so 
when  a short  duration  disturbance  occurs  In  the  Input  signal, 
the  effect  on  the  output  signal  will  die  out  after  a finite 
number  of  sampling  periods. 

The  frequency  response  of  a digital  filter  Is  periodic  and  can  be 
expanded  as  a Fourier  series  In  the  frequency  domain.  If  only  a finite 
number  of  terms  of  the  Fourier  series  Is  taken  frr  practical  purposes, 
we  must  truncate  the  Fourier  series.  Since  the  requency  response  of 
an  Ideal  filter  Is  always  specfflced  as  a discontinuous  function  with 
a point  of  discontinuity  occurring  at  the  cutoff  frequency,  the  trunca- 
tion of  Fourier  series  involves  the  Gibbs'  phenomenon  and  the  under- 
damped side  lobe  problems.  The  most  straightforward  way  to  solve  these 
problems  Is  through  the  use  of  windows. 

In  this  chapter,  we  will  develop  a concsc*'-  procedure  to  design 
nonrecurslve  filters  via  the  window  method^"^.  Starting  with  low- 
pass  filters  and  extending  the  discussion  to  high-pass,  band-pass  and 
band-rejection  filters,  a new  technique  of  frequency  transformation 
will  be  established  su:h  that  the  high-pass,  band-pass  or  band-rejection 
filters  designed  u«i»»g  this  technique  wilt  remain  nonrecurslve. 
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4.1  LOW-PASS  FILTER  DESIGN 


The  transfer  function  of  an  nth  order  nonrecursive  filter  can  be 


expressed  as 


F (Z)  « £ a z"k 

l._A  ^ 


Its  frequency  response  Is 


F(eJwT)  - £ a e'JkuT 


where  T Is  the  sample  period.  Assuming  n » 2m,  we  have 


F(.j«T)  . ,-J-T  £ ^ ,-Jk.T 


where 


bk  " ak-n 


Defining 


Then,  we  obtain 


2 £ bk  ,-ik»T 


F(eJ“T)  - e'JmuT  F (eJwT) 


In  fact,  Eq.  (5)  can  be  treated  as  the  Fourier  series  In  the  frequency 
domain  with  period  u>s  « which  is  the  sampling  frequency.  In  other 
words,  we  have 
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(3)  The  order  of  the  nonrecurslve  filter  Is  2m,  where  m Is  • 

(4)  Find  the  Inverse  Fourier  Transform  of  the  magnitude  response, 

| F (Z) I - |F  (Z)  I , by  using  the  FFT  algorithm. 

m 

(5)  The  first  m + 1 components  of  the  Inverse  Fourier  Transform 
of  |F(Z) | are  b^,  for  k ■ 0,  1,  ....  m. 

(6)  In  order  to  have  a smooth  frequency  response,  windows  must  be 
used  to  modify  the  coefficients  b^,  obtaining  b^,  for  k « 

0,  1,  • . , , m. 

(7)  The  transfer  function  of  the  nonrecurslve  filter  Is  given  by 

F(Z)  - £ a.  Z'k  (12) 

k-0  K 

where  a.  ■ b.'  _ and  b'.  ■ b!  . 

ic  k-m  k k 

t 

It  Is  often  necessary  to  find  the  transient  and  frequency  responses 
of  the  nonrecurslve  filter.  From  Eq.  (12),  the  impulse  response  is 

f(t)  - p ak  6 (t-kt)  (12a) 

Therefore,  the  Impulse  response  Is  very  easily  found  once  the  coefficients 
of  nonrecurslve  filters  are  given.  The  frequency  response  of  the  non- 
recurslve filter  Is  as  shown  In  Eq.  (2).  Of  course,  we  can  find  the 
frequency  response  directly  from  Eq.  (2)  If  the  coefficients  (a^'s)  ar* 
given.  However,  the  frequency  response  is  the  Fourier  Transform  of 
the  impulse  response  as  shown  In  Eq.  (12a).  Thus,  the  "fast"  way  to 
find  the  frequency  response  Is  to  take  the  Fourier  Transform  of  the 
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impulse  response  using  the  FFT  algorithm  (Appendix  A).  Program  No.  5 
(Appendix  C)  is  the  keyboard  program  for  nonrecursive  filter  design 
and  for  finding  the  frequency  response. 

Example  I.  A low-pass  filter  with  the  amplitude  response 


SI  | f | < 100  Hz 

0 | f | > 100  Hz 

and 

k - * -0.005  sec 

4>  du 

is  required.  The  sampling  frequency  is  1600  Hz.  Find  the  non recursive 
filter  using  the  window  method. 


Solution.  Since  k^  " “ -0.005  * -mT,  we  have 


m - - 8 


(13) 


Using  four  different  windows,  the  coefficients  of  four  different  non- 
recurslve  filters  were  determined  and  listed  in  Table  4- 1 . The 
frequency  response  of  each  of  the  resulting  filters  Is  shown  In 
Figures  4-1  through  4-4 . 

Figure  4-3  depicts  the  frequency  response  of  the  nonrecursive 
filter  designed  with  the  trapezoidal  window  (b  • 0.438);  note  that  the 
peak  of  the  highest  side  lobe  Is  -39  dB.  In  Figure  4-4  the  frequency 
response  of  the  nonrecursive  filter  designed  with  the  Hanning  window 
Is  shown;  note  that  for  this  filter  the  peak  of  the  highest  side  lobe 
Is  -44  dB.  From  Figures  4-3  and  4-4  we  see  that  the  trapezoidal  window 
and  the  Hanning  window  are  comparable;  however,  an  Important  advantage 
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FIGURE  1»-1.  FREQUENCY  RESPONSE  OF  LOW-PASS  FILTER  WITH  THE  RECTANGULAR  WINDOW 
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of  the  trapezoidal  window  is  that  it  is  much  easier  to  construct  than 
the  Hanning  window. 


TABLE  4-1 

COEFFICIENTS  OF  LOW-PASS  FILTER  DESIGNED  BY  FOUR  DIFFERENT  WINDOWS 


COEFFICIENTS 

WINDOWS 

RECTANGULAR 

TRAPEZOIDAL 

HANN 1 NG 

b = 0.375 

mm- 

a8 

0.1253 

0.1252 

0.1252 

0.1252 

a7  " a9 

0.1221 

0.1220 

0.1220 

0.1174 

a6  " a10 

0.1128 

0.1127 

0.0966 

0.0962 

a5  = a11 

0.0982 

0.0818 

0.0701 

0.0678 

a4  “ a12 

0.0796 

0.0531 

0.0455 

0.0398 

a3  " a13 

0.0588 

0.0294 

0.0252 

0.0182 

a2  = al4 

0.0374 

0.0125 

0.0107 

0.0055 

al  " a15 

0.0172 

0.0029 

0.0025 

0.0007 

a0  " al6 

-0.0003 

0 

0 

0 

***  4.2  FREQUENCY  TRANSFORMATION  I:  LOW-PASS  TO  LOW-PASS  TRANSFORMATION 

[18  23l 

Cons tant in  Ides 1 ’ developed  a set  of  formulas  for  the  frequency 

transformation  of  digital  filters.  Although  his  formulas  work  well  for 
recursive  filters,  they  are  not  suitable  for  nonrecurslve  filters.  When 
any  of  his  formulas  are  used  to  transform  a nonrecurslve  filter,  the 
resultant  is  no  longer  nonrecurslve. 

In  this  and  subsequent  sections,  we  will  develop  a set  of  trans- 
formations for  nonrecurslve  filters  to  construct  (1)  a low-pass  filter 
with  a different  cutoff  frequency,  (2)  a high-pass  filter,  (3)  a 
band-pass  filter,  and  (4}  a band-rejection  filter,  from  a low-pass  filter 
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prototype.  Let  us  start  with  the  low-pass  filter. 


Figure  4-5 (a)  shows  the  amplitude  response  for  the  ideal  low-pass 
filter  with  cutoff  frequency  f£  and  sampling  frequency  f . A new  non- 
recurslve  filter  with  cutoff  frequency  f^  as  shown  In  Figure  4-5 (b)  is 
required.  The  coefficients  of  the  low-pass  filter  with  a rectangular 
window  are  given  by 


b.  - — 
k »s 


juiTt  -jkuT 


F (eJUH ) e 


F (J2lTfT)  e’j2l,kfT  df 


where 


If  the  low-pass  nonrecurslve  filter  has  cutoff  frequency  f , the  coeffl- 

c 

dents  are 


. , ff«  . 


J2*kfT  df 


sin  (2kfcT) 


When  the  window,  w(t).  Is  used  for  smoothing  the  amplitude  response,  the 
modified  coefficients  become 


b.  ■ b.  * w. 
kc  kc  k 


(18) 


where  w^  Is  the  window  function  sampled  at  t « kT,  or  w^  = w(kT).  If 
the  cutoff  frequency  is  f^,  we  ha\ 
filter  with  a rectangular  window: 


the  cutoff  frequency  is  f^,  we  have  the  coefficients  of  the  nonrecursive 


. sin  (2-irkf  .T) 

■ — k-rr— 


(19) 


If  the  window,  w(t),  is  used  for  smoothing  the  amplitude  response,  the 
modified  coefficients  of  the  nonrecursive  filter  are 


bkd  " bkd  * wk 


(20) 


From  Eqs.  (18)  and  (20)  we  have 


b - kd  h 
bkd  . * kc 

bkc 

Substituting  Eqs.  (17)  and  (19)  into  Eq.  (21),  yields 


(21) 


sin  (2irkfdT) 
bkd  " sin  (2wkf  t)  bkc 


(22) 


If  sin  (2*kf  T)  i 0. 


Let  us  now  state  the  following  rule  for  the  low-pass  to  low-pass 
transformation: 


If  the  low-pass  nonrecursive  digital  filter  is  given  by 


F(Z)  - Z'm  £ b Z'k 
k— m kc 


(23) 


75 


f 


J 


with  cutoff  frequency  f , the  new  filter 


f(Z)  . £ bkd  z'k 


k—m 


where 


(24) 


sin  (2nkfdT) 
bkd  sin  (2-irkf cT)  bkc 


(25) 


will  have  the  cutoff  frequency  f^,  if  sin  (2nkfcT)  + 0. 

Example  2.  Using  the  results  of  Example  1,  design  a low-pass  nonrecursive 
filter  with  a cutoff  frequency  of  200  Hz  and  the  same  sampling  frequency 
(1600  Hz). 

Solution.  From  Eq.  (25),  we  have 


sin 


kd 


sin 


UT 


J 


(26) 


The  coefficients  of  the  required  nonrecursive  filter  are  listed  in 
Table  4-2.  The  frequency  responses  of  these  three  nonrecursive  filters 
are  shown  In  Figures  4-6,  4-7  and  4-8,  respectively. 

4.3  FREQUENCY  TRANSFORMATION  li:  LOW-PASS  TO  HIGH-PASS  TRANSFORMATION 


Figure  4-9  shows  the  amplitude  response  of  an  Ideal  high-pass  filter. 

f 

If  the  vertical  axis  Is  shifted  to  , we  find  that  the  amplitude 

response  of  the  high-pass  filter  with  cutoff  frequency  f^  and  sampling 

frequency  f becomes  the  amplitude  response  of  the  low-pass  filter  with 
s f 

cutoff  frequency  f ■ - f . 

c i n 
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FREQUENCY  RESPONSE  OF  LOW-PASS  FILTER  WITH  RECTANGULAR  WINDOW 


Mathematically  we  have 


Fh(eJ^ri)  = Fh  (2) 


Z •=  e 


j2irfT 


F£[,j2,T(f  - w 


where  Fh(Z)  Is  the  transfer  function  of  the  high-pass  filter  and  (Z) 
is  the  transfer  function  of  the  low-pass  filter.  Assume  that  the  trans 
fer  function  of  the  low-pass  filter  with  cutoff  frequency  f is 


From  Eq.  (27)  we  have 


F,{Z)  - £ , Z- 

k^m  K 


-Jirf  T, 
Fh(Z)  - F,(.  s Z) 


FIGURE  A-9.  THE  AMPLITUDE  RESPONSE  OF  AN  IDEAL  HIGH-PASS  FILTER 


Since  fs  « y • Eq.  (29)  becomes 


Fh(Z)  - F,(e'jlTZ) 
- F1 (-Z) 


Combining  Eqs.  (28)  and  (30),  yields 

Fh(Z)  - £ ak  (-Z) 

k— m 


(30) 


In  brief,  we  have  the  following  rule  for  the  low-pass  to  high-pass 
transformation: 


If  a high-pass  nonrecurslve  filter  with  — - k^,  cutoff  frequency 

f.  , and  sampling  frequency  f Is  required,  we  first  design  a low-pass 

S f 

nonrecurslve  filter  with  cutoff  frequency  f ■ tt-  - f.  and  4^-  = k •, 

c 2 h du>  4> 

then  substituting  Z by  -Z,  we  arrive  at  the  high-pass  nonrecursive  filter 
with  the  specified  performance. 

Example  3.  Design  a nonrecurslve  filter  with  a sampling  frequency  of 

1600  Hz , cutoff  frequency  of  400  Hz,  and  k ■ -0.005  seconds. 

$ 


Solution.  We  first  design  a prototype  low-pass  filter  with  cutoff  fre- 
fs 

quency  f£  ■ ^ f^  - 400  Hz,  k^  * -0.005,  and  m = ^ “ 8.  Using  the 

near  optimal  trapezoidal  window  with  b ■ 0.438,  we  determine  the  transfer 
function  of  the  prototype  low-pass  filter 


F,(Z)  - Z-8 


(32) 
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where 

bg  “ b.g  - 0 

b7  - b_7  - 0.0065 

b6  “ b-6  “ "O-000! 
b$  - b_5  - -0.0273 

bi,  * b-i,  “ 0.0002 

bj  - b_j  ” 0.0759 
b2  - b_2  - -0.0003 
b,  - b_,  - -0.3183 

bQ  - 0.5001  (35) 

The  amplitude  response  and  phase  characteristics  of  the  filter  defined 
by  Eq.  (34)  are  shown  in  Figure  4-10. 

4.4  DESIGN  OF  A NONRECURSIVE  BAND-PASS  FILTER 

Let 

Fh<»  ■ z”  £ v z'k  <*> 

and 

Ft(Z)  " z’m  bk  Z"k  (37) 

be  the  transfer  functions  of  a high-pass  filter  with  cutoff  frequency 
f^  and  a low-pass  filter  with  cutoff  frequency  fc>  respectively. 
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PHASE 


Figure  1 1 (a)  shows  the  amplitude  responses  of  Eqs.  (36)  and  (37). 
Assume  that 


f > 
c 


Then,  the  amplitude  response  of 


(38) 


Fb(Z)  - Fh(Z)  Ft(Z)  (39) 

as  shown  In  Figure  4- 1 1 (b) , exhibits  the  band-pass  characteristic,  with 
pass  band  f j <_  f <_  f2>  where  fj  - fh  and  - f£.  Therefore,  we  can 
construct  the  band-pass  nonrecurslve  filter  by  cascading  a low-pass  and 
a high-pass  filter  with  proper  cutoff  frequencies.  Substituting  Eqs.  (36) 
and  (37)  into  Eq.  (39) » yields 

F‘,z)  ■ z'2l£  b*  z1  ,w 


Eq.  (40)  Indicates  that  the  band-pass  nonrecurslve  filter  has  a linear 
phase  relationship  with  - -2mT. 

Now,  we  can  state  the  following  rule  for  the  construction  of  a band- 
pass filter: 


and  F*2  (Z) , with  cutoff  frequencies  ip  - fj  end  f2,  respectively,  where 

<*♦.  d*2  k* 

f#  Is  the  sampling  frequency;  and  the  phase  slopes  gp- ■ gp-  • j3-  , 
where  and  $2  are  the  phase  functions  of  (2)  and  Fjt2 (Z) , then  the 


transfer  function  of  the  desired  band-pass  filter  Is  given  by 


Fb(Z)  - Ft,(-Z)  F£2(Z) 


(41) 
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Example  A.  Design  a band-pass  filter  with  pass  band  200  Hz  < f < 1*00  Hz 


and  ■ -0.01  sec.  The  sampling  frequency  is  1600  Hz. 


Solution.  Since  f,  » 200  Hz,  f,  - 1*00  Hz  and  f ■ 1600  Hz,  we  first 

1 s f 

design  two  low-pass  filters  with  cutoff  frequencies  ~ - f j ■ 6°°  Hz 

and  400  Hz,  respectively,  and  k,.  ■ k.«  ■ -0.005.  The  order  of  low-pass 
i,  9*  9* 

'e  1 

filters  Is  in  ■ ■ 8.  Following  the  procedures  given  In  Section  l*.l 

and  using  the  near  optimal  trapezoidal  window  with  b - 0.1*38,  yields 


F,,  (?)  - 2- 


FtJ(2)  - 2- 


where 


a8  " a-8 


a7  * *-7 


*6  " *-6 


*5  * *-5 


•l.  * *-l* 


*3  " a-3 


a2  " *-2 


al  " *-l 


-0.005 


0.016 


-0.020 


-0.001 


0.051* 


-0.137 


0.225 


0.750 


l 


and 


by  - b_7  - -0.0065 

bg  - b_g  ■ -0.0001 

b5  " b_5  - 0.0273 
b^  ■ b-l4  ■ 0.0002 
b3  - b_3  - -0.0759 
b2  - b_2  - -0.0003 
b,  - b_,  - 0.3183 

bQ  - 0.5001  (45) 

Now,  using  Eq.  (41)  the  transfer  function  of  the  band-pass  filter 
becomes 

Fb(2)  * Ftl(_Z)  Ft2(2) 
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where 


- 0 
- 0 

- -0.0001 
- -0.0002 
- -0.0001 

- 0.0004 

- 0.0005 

- -0.0004 

- -0.0001 

- -0.0015 

- 0.0145 

- 0.0460 

- -0.0001 

- -0.1320 

- -0.1429 

- 0.0935 

- 0.2410 


(47) 
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The  frequency  response  of  the  band-pass  filter  as  expressed  In  Eq.  (46) 
Is  shown  In  Figure  4-12. 


4.5  DESIGN  OF  A NONRECURSIVE  BAND-REJECTION  FILTER 
Let 

Fh(z>  ■ z'"  H z'k  <w> 

and 

F£(Z)  - Z'm  bk  Z'k  (49) 

be  the  transfer  functions  of  a high-pass  filter  with  cutoff  frequency 

f.  and  a low-pass  filter  with  cutoff  frequency  f , respectively.  Figure 
n c 

4-1 3 (a)  shows  the  amplitude  responses  of  Eqs.  (48)  and  (49).  Assume 
that 

f.  > f (50) 

n c 

Then  the  amplitude  response  of 

Fr(Z)  - Fh(Z)  + FA(Z)  (51) 

as  shown  In  Figure  4-13 (b) . exhibits  the  band- reject  Ion  characteristics, 
with  rejection  band  f j < f < f^  where  f j - ffi  and  f2  - f^.  Therefore, 
we  can  construct  the  band- reject I on  nonrecurslve  filter  by  the  parallel 
connection  of  a high-pass  and  a low-pass  filter  with  proper  cutoff  fre- 
quencies. Substituting  Eqs.  (48)  and  (49)  Into  Eq.  (51),  yields 

F (Z)  - I-  (ak  ♦ bk)  z'k  (52) 

r l£Pm 
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BAND-PASS  FILTER  WITH  TRAPEZOIDAL  WINDOW,  b - O.A38 


I 


f 

s 


(b)  BAND-REJECTION  FILTER 

FIGURE  4-13.  CONSTRUCTION  OF  THE  BAND-REJECTION  FILTER  FROM  A LOW-PASS 

AND  A HIGH-PASS  FILTER 
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It  Is  evident  that  the  band-rejection  filter  as  represented  by  Eq.  (52) 
has  a linear  phase  relationship  with  ^ ■ -mT. 

Accordingly,  we  can  now  state  the  following  rule  for  the  construction 
of  the  band- reject  Ion  filter. 

If  the  band- reject  I on  nonrecurslve  filter  with  rejection  band 
fj  if  i fj  and  ^ “ k^  is  required,  we  first  design  two  low-pass  filters 

(Z)  and  F42(z)  w*th  cutoFF  frequencies  fj  and  - f2,  respectively, 

d$. 

where  f is  the  sampling  frequency;  and  the  phase  slopes  are  -r — * 
d<|>2  5 dw 

■ k^,  where  and  $2  are  phase  functions  of  F4)(z)  and  F^Z),  then 

the  transfer  function  of  the  desired  band-rejection  filter  is  given  by 

Ff(Z)  - F41(Z)  + Ft2(-Z)  (53) 


Example  5.  Design  a band-rejection  filter  with  rejection  band  100  Hz  <_ 
f <^600  Hz  and  " -0.005  seconds.  The  sampling  frequency  is  1 600  Hz. 


Solution.  Since  fj  - 100  Hz,  f2  ■ 600  Hz  and  f - 1600  Hz,  we  first 

design  two  low  pass  filters  with  cutoff  frequencies  of  100  Hz  and  of 
f 

- f2  ■ 200  Hz,  respectively,  and  k^  - k^2  ■ -0.005  seconds.  There- 
fore, the  order  of  the  low-pass  filters  is  m ■ 8.  If  the  near  optimal 
trapezoidal  window  with  b - 0.438  Is  used,  we  have 


and 


F*>(z,-z'8  kfc  **  z‘k 

F«2|Z>  * Z'8  kte  bk  2'k 


(54) 


(55) 
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f 


where 


*8  " a-8 

0 

*7  “ a-7 

0.0025 

*6  “ a-6 

0.0107 

*5  * a-5 

0.0252 

a*  - •.<, 

0.0*55 

*3  “ *-3 

0.0701 

*2  * a-2 

0.0966 

a1  “ a-l 

0.1220 

ao 

0.1252 

(56) 

b8  " b-8 

0 

b7  " b-7 

0.00*6 

NO 

1 

XX 

» 

NO 

XX 

-0.0152 

b5  - b-5  ■ 

0.019* 

bl,  * b.i,  * 

0 

b3  * b-3  ' 

-0.0535 

b2  " b-2  ' 

0.136* 

b1  " b-1  ' 

-0.2253 

bn  1 

0.2502 

(57) 
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Therefore,  the  required  band-rejection  nonrecurstve  filter  is  given  by 


where 


Fr(z“’) 


(58) 


c8  " c-8  “ 0 

Cy  * C_J  ■ 0.0070  C. 

Cg  ■ c_g  » -0.0045 
Cj  » c_j  » 0.0446 
” c_^  » 0.045^ 

Cj  » c_^  *■  0.0167 
c2  * c-2  “ 0.2.129 

c1  - c_j  * -0.1032 

cQ  - 0.3754 

The  frequency  response  of  this  band-rejection  nonrecurslve  filter  5s 
shown  In  Figure  4-14. 

4.6  REMARKS 

The  design  of  various  kinds  of  nonrecurslve  digital  filters  using 
the  FFT  algorithm  and  windows  has  been  presented.  It  Is  evident  that  an 
arbitrary  frequency  response  can  be  realized  through  a nonrecurslve 
digital  filter  of  the  proper  order.  This  technique  will  be  employed  to 
design  nonrecurstve  digital  compensators  for  control  systems. 
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CHAPTER  V 

ANALYSIS  OF  DISCRETE-TIME  CONTROL  SYSTEMS* 

5.0  INTRODUCTION 

To  understand  and  analyze  a discrete-time  control  system,  we  often 
manipulate  and  transform  the  system  block  diagram,  and  then  study  the 
system  response  and  stability. 

In  this  chapter,  we  will  develop  a consistent  set  of  rules  for  dealing 
with  system  block  diagrams.  Block  diagram  operations,  simplifications, 
and  closed-loop  transfer  function  evaluations  will  be  derived  from 
direct  applications  of  these  rules;  i.e.,  the  Ragazzlnl-Zadeh  identity. 

As  far  as  the  transient  response  and  the  frequency  response  are 
concerned,  algorithms  have  been  developed  which  are  Fast  Fourier  Trans- 
form oriented;  therefore,  the  techniques  devised  are  to  be  used  with  the 
aid  of  minicomputers  or  microcomputers. 

Finally,  a novel  stability  criterion  based  on  the  return  difference 
Is  derived,  explained,  and  demonstrated  through  various  examples.  This 
criterion  Is  also  of  the  FFT-orlented  and  minicomputer-aided  type. 

In  general,  this  chapter  presents  methods  for  quantitative  and 
qualitative  evaluation  of  the  performance  of  a discrete-time  control  sys- 
tem. Based  on  a thorough  understanding  of  the  principles  presented  in 
this  chapter,  we  will  establish  the  procedures  for  the  design  work 
presented  in  the  following  chapter. 

5.1  FREQUENCY  RESPONSE  OF  DISCRETE-TIME  CONTROL  SYSTEMS 

If  the  transfer  function  of  a discrete-time  control  system  Is  given 
by 
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r 


an  + a.  z + a,  zZ  + ...  + a z 
G(z)  - -2 1 2 m 

bQ  + b.  z + b,  zZ  + ...  + b zn 


m 


(1) 


then  the  frequency  response  of  the  system  is  given  by 


G(jw)  - G(z) 


-eJwT 


eQ  ♦ a,  e^T  ♦ a,  e^“T  e^T 

b0  + b,  eJwT  + b2  ej2wT  + ...  ♦ bR  ejnwT 


(2) 


where  T Is  the  sampling  period  and  u «■  ~L  | s the  sampling  frequency. 
j(u>+u>s)T  * T 

Since  e • eJ“  , we  see  that  G(ju)  is  a complex  periodic  function 

In  the  frequency  domain  with  period  ws.  Therefore,  Eq.  (12)  can  be 
expressed  as  a Fourier  series; 


w 

G(jo»)  - ^ an  * 


■jnwT 


(3) 


n*-» 


The  Inverse  Fourier  Transform  of  G(ju),  which  Is  the  Impulse  response  of 
the  discrete-time  control  system.  Is  given  by 

g(t)  - F_1 [G(jw) J 


m 

- 23  «n  S(t-nT) 


n*-« 


<*> 


For  a causal  system,  the  Impulse  response  must  be 


g(t)  • 0 for  t < 0 


(5) 
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In  other  words,  In  Eq.  (A) 


a_  m 0 for  n < 0 
n 


For  simplicity,  Eqs.  (3)  and  (A)  can  be  rewritten  as 


and 


(6) 


(7) 


(8) 


The  conventional  methods  to  obtain  the  an  coefficients  of  Eq.  (8)  are 
the  partial  fraction  expansion,  long  division  algorithm,  etc.  A new 
method  of  using  the  OFT  to  find  the  frequency  response  G(jw)  and  then 
taking  the  inverse  Fourier  Transformation  of  G(jw)  to  get  g(t)  has  been 
developed  and  can  be  adopted. 

in  the  preceding  chapter,  the  frequency  response  of  the  nonrecurslve 
filter  was  determined  using  the  FFT  algorithm.  Eq.  (1)  can  be  rewritten 

as 

'<«>  ■ rel  <»> 

where 

N(z)  - am  zm"n  ♦ am_j  2m"n-l  + ...  + z'"n  + a0  z"n  (iO) 
and 

D(i)  ■ bfl  ♦ bn-)  z '+...+  b(  z*  n + bfl  i " (11) 
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and  m < n.  Of  course,  N(z)  and  D(z)  can  be  treated  as  transfer  functions 
of  two  nonrecurslve  filters.  Therefore,  utilizing  the  technique  developed 
In  Chapter  IV,  the  frequency  response  of  N (z)  and  D(z)  can  be  determined 
using  the  FFT  algorithm.  Next,  dividing  N(e^uT)  by  D(e-*wT),  the  frequency 
response  of  the  system  transfer  function  G(z)  Is  obtained.  The  procedure 
suggested  here  Is  actually  a fast  method  to  calculate  the  frequency 
response  of  discrete-time  control  systems.  Program  No.  6 (Appendix  C) 
was  developed  according  to  this  procedure. 

Example  I . Find  the  frequency  response  of  the  system  with  the  transfer 
function 


G(z)  - 


1 . 26Az 

0. 104z  + 0.368 


and  sampling  period  T ■ 0.1  seconds. 


Solution.  Rewrite  Eq.  (12)  as 


G(z)  • 


1 .264z 


1 - 0. 104z"  + 0.368z' 


Case  (1).  Choose  N - 64  for  performing  the  FFT. 

The  frequency  resolution  of  the  frequency  response  is  given  by 

tMmJTmW  (u) 

In  this  case,  T ■ 0.1  seconds  and  Aw  ■ jg-  ■ 0.9818.  The  resulting  frequency 

“s 

response  Is  shown  In  Figure  5-1  for  0 £ « £ • 


MUM 


Case  (2).  Choose  N - 1024  for  performing  the  FFT. 

In  this  case.  Aw  ■ ffr  “ 0-06136  Is  much  smaller  than  that  In 

w 

Case  (1).  The  frequency  response  for  0 <_  w < Is  shown  In  Figure  5-2. 
Increasing  the  value  of  N results  In  a much  smoother  curve  compared  to 
that  of  Case  (1). 

5-2  TRANSIENT  RESPONSE  OF  DISCRETE-TIME  CONTROL  SYSTEMS 

If  the  discrete-time  control  system  has  the  transfer  function  G(z) 
and  input  R(z),  then  the  output  is  given  by 

C(z)  - G(z)  R(z)  (15) 

and  the  transient  response  is 

c(t)  - Z-1[G(z)  R(z)]  (16) 

Many  methods  exist  for  finding  the  transient  response  c(t),  such  as  the 
long  division  method,  and  partial  fraction  expansion,  etc.  However,  If 
we  take  advantage  of  the  FFT  algorithm  for  finding  the  inverse  Z-transform 
as  In  Eq.  (16)  and  using  Eo.  (15),  we  have  the  Fourier  Transform  of  C(z)  \ 

which  can  be  expressed  as 

C(eJuT)  - C(eJwT)  R(eJuT)  (17) 

Then,  the  transient  response  Is 

c(t)  - F’1[C(eJwT)] 

- F^tCfe^1)  R(eJwT)]  (18) 
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Thus,  we  first  find  the  frequency  response  G(e^uT)  and  the  Fourier  Trans- 
form R(e^“T)  using  the  FFT,  as  discussed  In  Section  5.1,  then  we  use  the 
FFT  to  evaluate  the  inverse  Fourier  Transform  of  Gte-*^)  R(eJU,T)  to  deter 
mine  c(t).  Details  of  this  method  are  given  In  Appendix  B. 

5.2.1  UNIT  IMPULSE  RESPONSE 

The  unit  Impulse  input  Is  defined  as 

R(z)  - Z[6 (t) ] - 1 (19) 

From  Eq.  ( 1 8) , we  have  the  unit  Impulse  response 

c ( t)  - F',[G(eJuT)]  (20) 

Program  No.  7 (Appendix  C)  was  written  to  evaluate  the  unit  impulse 
response  from  the  transfer  function  G(z). 

Example  2.  Find  the  unit  Impulse  response  of  the  system  with  the  trans- 
fer function 

G(z)  ■ —x (21) 

z - 0.l04z  + 0.368 

and  sampling  period  T « 0.1  seconds. 

Solution.  Choose  N ■ 1024.  Utilizing  Program  No.  7 (Appendix  C)  the 
unit  impulse  response  is 
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c (t)  - 1.26A  6 (t-0. 1)  + 0.132  6(t-0.2) 

- 0.0452  6 (t-0. 3)  - 0.096  6(t-0.4) 
+ 0.157  6(t-0.5)  + 0.052  6 (t-0.6) 

- 0.053  6(t-0.7)  - 0.025  «(t-0.8) 

+ 0.017  5(t-0.9)  - 0.011  «(t-1.0) 

- 0.005  6 ( t— 1 . 1) 

The  resulting  transient  response  is  shown  In  Figure  5~3. 


FIGURE  5-3.  THE  UNIT  IMPULSE  RESPONSE  OF  EXAMPLE  2 


If  the  system  has  poles  at  z ■ 1 , the  frequency  response  ts  undefined 
In  the  neighborhood  of  u - 0;  therefore,  the  Impulse  response  cannot  be 
evaluated  directly  by  FFT  algorithm.  To  overcome  this  difficulty,  we 


can  factor  out 


and  define 


(23) 


Mz)  - (iT-L)k  G(z> 

such  that  no  pole  exists  at  z - I In  G,(z).  The  impulse  response  of 
G,(z)  is  determined;  and  the  method  for  evaluating  the  step  response, 
which  is  discussed  In  the  following  section.  Is  used  to  find  the  impulse 
response  of  G(z). 

5.2.2  UNIT  STEP  RESPONSE 


When  the  Input  is  a unit  step  function,  defined  as 


R(z) 


(24) 


Eq.  (24)  can  be  expanded  as 

R(z)  "l+z*+z2+...+zn+...  (25) 

Thus,  the  unit  step  response  is 
C(z)  ■ R(z)  G(z) 

- G(z)  + z_1  G(z)  ♦ z'2  G(z)  ♦ ...  + z‘n  G(z)  + ...  (26) 


and  we  have 


where 


c(t) 


g(t-iT) 


g(t)  - z",[G(z)] 


(27) 


(28) 
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Is  the  unit  impulse  response  of  the  system.  Since  the  physical  system 
Is  causal , 


g(t)  • 0 for  t < 0 


(29) 


From  Eq.  (27) • we  find 

n 

c(nT)  - 2 9(nT)  for  n ■ 0,  1,  ...  (30) 

1-0 

The  Fourier  Analyzer  has  an  Integration  algorithm  which  actually  performs 
the  lummation  exactly  as  shown  In  Eq.  (30);  therefore,  the  Integration 
operation  can  be  used  to  evaluate  the  unit  step  response.  Program  No.  8 
(Appendix  C)  was  written  to  determine  the  unit  step  response  from  the 
transfer  function  of  the  system. 


Example  3.  Find  the  unit  step  response  of  the  system  with  transfer 
function 

/ 

GW--- 1^2 (31) 

z - 0.104z  + 0.368 

and  sampling  period  T - 0.1  seconds. 

Solution.  Using  Program  No.  8 (Appendix  C),  the  unit  step  response  is 


c (t) 


- 1 .264  6(t-0.1) 

+ 0.941*  6(t-0.3) 
+ 1.005  «(t-0.5) 
+ 1.004  6(t-0.7) 
+ 1.007  «(t-0.9) 
+ 0.998  6( t-l.l) 


• 1.395  «(t-0.2) 

+ 0.849  6 (t-0.4) 

+ 1.056  6 (t-0.6) 
+ 0.980  «(t-0.8) 

♦ 1.002  5(t-1.0) 
+ ... 


(32) 
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FIGURE  5 -k.  THE  UNIT  STEP  RESPONSE  OF  EXAMPLE  3 


Example  4.  Find  the  unit  Impulse  and  unit  step  responses  of  the  system 
with  the  transfer  function 

S(.) '•?!.  (33) 

(z  - l)(z  - 0.5Z  + 0.5) 

and  sampling  period  T - 0.1  seconds. 


Solution.  First,  define 


G,(z) 


z - I 1.5  (z  ♦ I) 

2 zi  - 0.5z  + 0.5 


The  unit  Impulse  response  of  G,(z)  Is 


(34) 


>8 


g,(t)  - 1.501  6 (t-0.  1)  + 2.25  6 (t-0.2) 

+ 0.0375  <5(t-0.3)  - 0.938  6(t-0.4) 

- 0.657  6 (t-0. 5)  + 0.141  6 (t-0. 6) 

+ 0.399  «(t-0.7)  ♦ 0.129  6 (t-0.8) 

+ 0.  135  6(t-0.9)  + 0.002  Mt-1.0) 

+ 0.067  6 (t-1 . 1 ) + ...  (35) 

From  Eq.  (30).  the  unit  Impulse  response  of  G(z)  Is 

g(t)  - 1.501  fi(t-O.l)  + 3.751  «(t-0.2) 

+ 4.126  6(t-0.3)  + 3.188  6(t-0.4) 

♦ 2.532  «(t-0.5)  + 2.673  6(t-0.6) 

+ 3.071  «(t-0.7)  ♦ 3.200  6(t-0.8) 

+ 3.065  6(t-0.9)  + 2.933  6(t-1.0) 

+ 2.935  «(t-1.1)  + ...  (36) 

and  the  unit  step  response  of  G(z)  Is 

h(t)  - 1.51  6 (t-0. 1 ) + 5.26  6(t-0.2) 

+ 9-38  6(t-0.3)  + 12.57  «(t-0.4) 

+ 15.10  6(t-0.5)  + 17.77  6 (t-0. 6) 

+ 20.85  fi(t-0.7)  + 24.05  6(t-0.8) 

+ 27.11  6 (t-0.9)  + 30.04  6 (t-1 .0) 

+ 32.98  6(t-1.1)  + ...  (37) 

Both  g(t)  and  h(t)  are  shown  In  Figures  5“5  and  5“6  respectively. 


Example  5.  Find  the  frequency,  unit  Impulse,  and  unit  step  responses 
of  the  discrete-time  system  with  the  transfer  function 


where 


aQ  - -0.006 

bQ  - 0.006 

aj  ■ 0.044 

bj  - -0.044 

a2  - -0.262 

b2  - 0.262 

a^  - -0.516 

b3  - 0.516 

- 0.152 

b^  - -0.152 

a?  - -0.210 

b$  - 0.210 

a6  - 1.116 

b6  - -1.116 

a7  - 1.682 

b?  - -1.682 

b8  " *-° 

and  sampling  period  T - 0.5  seconds 

• 

Solution.  Utilize  Program  No.  8 (Appendix  C)  to  find  both  frequency  and 

transient  responses  of  the  system. 

Figure  5"7  shows  the 

resulting  fre- 

quency  response. 

The  unit  impulse 

response  is  given  by 

g(t) 

- 0.4205  6 (t— 0.5) 

f 0.4458  6 (t-1 .0) 

+ 0.2566  6 (t-1 .5) 

+ 0.2510  6(t-2.0) 

+ 0.0403  6 (t-2.5) 

- 0.0290  6 (t-3.0) 

- 0.0797  «(t-3.5) 

- 0.0939  6(t-0.4) 

- 0.1035  «(t-4.5) 

- 0.0861  «(t-5.0) 

- 0.0597  «(t-5.5) 

- 0.0350  fi(t-6.0) 

- 0.0139  6 (t-6.5) 

+ 0.0030  6(t-7.0) 

+ 0.0139  «(t-7.5) 

♦ 0.0185  6 ( t-8 . 0) 

+ 0.0186  6 (t-8.5) 

+ 0.0159  6 ( t— 9 • 0) 

+ 0.0117  «(t-9.5) 

+ 0.0070  «(t-10.0)  + 

(39) 

m 


FIGURE  5-9. 


and  the  unit  step  response  is 

Mt)  - 0.421  «(t-0.5)  ♦ 0.876  6 (t- 1 .0) 

♦ 1.133  «<t-1.5)  + 1.384  6 (t-2.0) 

♦ 1.424  «(t-2.5)  ♦ 1.395  «(t-3-0) 

♦ 1.315  6(t-3-5)  + 1.222  «(t-4.0) 

+ 1.118  $(t-4.5)  ♦ 1.032  «(t-5.0) 

+ 0.973  «(t-5.5)  + 0.939  «(t-6.0) 

+ 0.924  5 (t-6.5)  + 0.927  «(t-7.0) 

♦ 0.941  «(t-7.5)  + 0.959  «(t-8.0) 

+ 0.978  6(t-8.5)  + 0.994  <$(t-9.0) 

+ 1.005  6(t-9.5)  + 1.012  6(t-10.0)  ♦ ...  (40) 

Both  g(t)  and  h(t)  are  shown  In  Figures  5-8  and  5-9,  respectively. 

5.3  BLOCK  DIAGRAM  OPERATIONS  FOR  DISCRETE -TIME  CONTROL  SYSTEMS 

To  find  the  closed  loop  transfer  function  of  a discrete-time  system, 
we  must  first  derive  a set  of  block  diagram  operations  for  simplifying 
the  system  structure.  The  following  outlines  some  elements  and  funda- 
mental operations  for  simplifying  the  block  diagram  of  a discrete-time 
system: 

e The  S -domain 


(1)  Components 


(a)  CC1  - Dynamic  device 


X(s)  - G(s) 


Y (s)  - G(s)  X (s) 


(b)  CC2  - Summing  device. 

X,(s) Y(»)  . x,(s)  ± X2(s)  ± ...  ± Xn(s) 


x2(s)' 


X (s) 

n 


• The  Z-domaln. 


(1)  Components. 

(e)  DC1  - Dynamic  device. 


(c)  D03  - Feedback  operation. 


e The  hybrid  domain. 
(I)  Component. 


(a)  HC 1 - Sampler  (discrete  input). 

X(z) -(T) * x(z) 

(b)  HC1  - Sampler  (continuous  Input). 

X(s) -/T) ►Y(z)  - Z(X(s)] 


(2)  Operation. 


HOI  - Ragazzlni-Zadeh  identlty^'^. 

X(s) G ( s ) — ^Y(z)  - Z[G  (s)  Z [X  (s)  ] ] 

W _TV_y  - Z[G(s)]  Z[X(s)] 

All  of  the  previous  operations  are  based  on  the  linear  graph  theory  with 
the  exception  of  the  Ragazzinl-Zadeh  Identity.  Through  the  use  of 
these  five  components  and  seven  fundamental  operations,  we  can  derive 
other  useful  operations: 

• Insert  or  delete  the  sampler  - The  sampler  can  be  Inserted 
Into  or  deleted  from  the  path  through  which  the  digital  sig- 
nal Is  passed.  This  is  based  on  the  property  of  the  sampler 
(HCI). 


This  operation  results  directly  from  Ragazzlni-Zadeh  identity. 
• Moving  a summing  device  beyond/ahead  of  devices  or  samplers. 
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Moving  the  take-off  point  beyond/ahead  of  devices  or  samplers. 


5.4  CLOSED  LOOP  TRANSFER  FUNCTION 

The  block  diagram  operations  given  In  the  previous  section  will  be 
used  to  determine  the  closed  loop  transfer  function  of  systems  presented 
In  this  section.  Consider  the  system  shown  In  Figure  5“ 10(a) ; moving  the 
summing  device  beyond  the  first  samplier  and  the  take-off  point  ahead 
of  G2(s) , the  block  diagram  transforms  into  that  shown  In  Figure  5~ 1 0 (b) . 
Next,  using  the  Ragazzinl-Zadeh  Identity  (Operation  HOI)  to  eliminate  the 
samplers,  we  have  Figure  5~ 1 0 (c) . Finally,  applying  the  feedback  opera- 
tion 003  and  the  series  operation  D01,  the  transfer  function  can  be 
expressed  as 


F(z)  ’ £{fi‘  rr 


Gj(z)  G 2(z) 


where 


G2H(z)  - z|g2(s)  H(s)] 


Consider  the  missile  launching  system  with  a digital  controller  as 
shown  in  Figure  5-1 1(a).  The  launcher  dynamics  are  given  by 


ci(,)  ' 7TTV 


The  zero-order  hold  Is 


Go(s)  " ■ V-5  m 

where  T Is  the  sampling  period  (T  • 0.1  seconds).  The  digital  compensator 
is  given  by 

D (2)  - 0-0577(z  - 0.92) (z  - 0.716)  (*5) 

c z2  + 0.22z  + 0.41*6 

Moving  the  take-off  point  beyond  the  last  sampler  In  the  feed-forward 
path  and  eliminating  the  first  sampler,  the  block  diagram  reduces  to 
Figure  5-11  (b).  Applying  the  Ragazzlni-Zadeh  identity  (Operation  HOI) 
yields  Figure  5~11(a)  where 

G(z)  - z[gq(s)  Gl(s)]  (46) 

Since  Gq(s)  is  given  in  Eq.  (44),  we  can  express  G(z)  as 

0(r)  - z[icL(,)]  - z[isL(.)  «'Ts] 

- 0 - i)  z[iy»]  M) 

In  other  words,  we  can  write  the  zero-order  hold  as 

G0(s,z)  - 1^-  (48) 

If  we  replace  GQ(s)  In  Eq.  (46)  by  Eq.  (48)  and  then  apply  the  Ragazzlni- 
Zadeh  Identity,  we  can  Immediately  obtain  the  result  of  Eq.  (47).  From 
Eq.  (47)  we  have 
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(*9) 


G(z) 


_ 3.925  x lo'^z  + 3.35)  (z  + 0.26) 
* fz  - l)(z  -■  0.901*8)  (z  - 0.8187) 


Therefore,  the  closed  loop  transfer  function  is  given  by 
Dc(z)  G(z) 

F(z)  " \ ♦ 0c( z)G{z) 


. 3.925  x lO'Nz  + 3.35)  (z  + 0.26)  (z  - 0.92)  (0.0577z  - 0.403) 
z2(z  - 0.93*0  (z2  - 1.508z  + 0.6065) 


(50) 


5.5  STABILITY  OF  DISCRETE-TIME  CONTROL  SYSTEM 

Numerous  criteria,  such  as  the  Schur-Cohn  criterion,  the  Inner 
theory,  the  Ronth-Hurwi tz  criterion,  the  Nyquist  criterion,  etc.,  are 
available  for  testing  the  stability  of  discrete-time  control  systems. 

Most  of  these  criteria,  with  the  exception  of  the  Nyquist  criterion, 
are  algebraic  algorithms  and  are  based  on  the  closed  loop  characteristic 
equation.  The  original  Nyquist  criterion^  is  based  on  the  encircle- 
ment of  the  frequency  plot  of  the  return  ratio  function  of  the  closed 
loop  system.  Since  the  frequency  response  of  the  discrete-time  control 
system  can  be  found  utilizing  the  FFT  algorithm,  It  is  evident  that  the 
Nyquist  criterion  can  be  an  effective  tool  for  stability  testing  via 
the  FFT. 

[281 

Recently,  a generalized  Nyquist  criterion1  1 based  on  the  return 
difference  rather  than  the  return  ratio  has  been  derived  and  successfully 
used  for  testing  the  stability  of  lumped  and  distributed  parameter  con- 
tinuous systems.  Subsequently,  we  will  derive  Its  counterpart  for 
discrete-time  systems. 

If  a discrete-time  control  system  has  a feed-forward  path  transfer 
function  G(z)  and  a feedback  path  transfer  function  H(z),  the  closed 
loop  system  transfer  function  Is  given  by 
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c/_\  G(z) 

F'  ' 1 + G(z)  H (z) 


(51) 


Similar  to  the  definitions  for  continuous  system,  we  can  define  the 
return  ratio  of  the  discrete-time  system  as 


Q(z)  - G(z)  H(z) 


(52) 


and  the  return  difference  as 


P(z)  - 1 + Q(z) 


(53) 


Let  us  express  the  return  ratio  as 


Q(z) 


A(z) 

bTzT 


(5M 


where  A(z)  and  B(z)  are  relative  prime  polynomials,  and  the  order  of 
A(z)  Is  not  greater  than  that  of  B(z).  Substituting  Eq.  (51*)  into  Eq. 
(51).  yields 


G(z)  B(z) 

F(z)  “ A(z)V  B(zT 


(55) 


Assuming  that  there  Is  no  common  factor  between  the  denominator  of  G(z) 
and  the  numerator  of  H(z),  then  G(z)  B(z)  In  Eq.  (55)  Is  a polynomial; 
B(z)  Is  the  open  loop  characteristic  function  and  A(z)  + B(z)  Is  the 
closed  loop  characteristic  function  Aj(z).  Thus  the  return  difference 
becomes 

..v  A(z)  + B(z)  _ A1 (z) 

p(z)  - 6{zl° — ettt 


closed  loop  characteristic  function 
open  loop  characteristic  function 


(56) 
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The  original  Nyquist  criterion  for  discrete-t Ime  systems  Is  to  count  the 

number  of  encirclements  of  the  frequency  response  of  the  return  ratio 
w u 

fron,  _L  to  - •—  , where  u>s  is  the  sampling  frequency,  about  the  critical 

point  (-1  + jO).  The  generalized  Nyquist  criterion  considered  here  Is 

to  count  the  total  phase  change  of  the  frequency  response  of  the  return 
w 

difference  from  0 to  ■—  ; therefore,  the  new  criterion  is  much  simpler 
than  the  original  one. 


Before  preceding  any  further,  the  following  definitions  are  required: 

• Definition  1.  The  open  arc  on  the  unit  circle  with  the 

center  at  the  origin  of  Z-plane,  from  Uj  to  u>j+j  is  the 

ju,T  Jwj+iT 

open  segment  of  arc  (e  , e ),  where  T Is  the 

sampling  period. 


Definition  2.  The  phase  angle  change  of  the  mapping  of 
Jw.T)  + 

an  open  arc  l(l)  ■ (e  , e ) by  a certain  complex 

function  P(z)  Is  defined  by 


A6 


*0) 


(v, )•(»:) 


(57) 


where 


;♦<  ■ ’ c))\ 


(58) 


and 


{*"'  P(,J<"'T  * °)} 


(59) 
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• Definition  3-  The  total  phase  change  of  the  mapping  of 
all  open  arcs  on  the  upper  half  of  the  unit  circle  by 
P(z)  Is  defined  as 

ae  - £ *et(1)  m 

where  Is  the  phase  angle  change  due  to  the 

mapping  of  the  I th  open  arc. 

Assume  that  the  return  difference  P(z)  has  the  following 
properties: 

(1)  P(e^w^)  4 0 for  0 <_  id  <_  . 

(2)  There  are  6 poles  of  P(z)  on  the  unit  circle  and 
Y poles  out  of  the  unit  circle,  where  0 and  y 
are  finite  Integers. 

(3)  P*(z)  * P(z*)  where  * means  complex  conjugate. 

Consider  the  Nyqulst  contour  for  discrete-time  system  as  shown  In 
Figure  5-12.  From  the  principle  of  argument,  It  Is  known  that  the  total 
phase  change  of  P(z)  along  the  contour 

r - l(-n)  Uc(-n+1)  ...  Ut(-l)  Uc(0)  U*(1)  U ... 

...  U£(n)  Uc(n)  (61) 

Is  given  by 

Aey  - 2n[(NA-M)  - (Nr-  y)]  (62) 
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where  N^  Is  the  order  of  Aj(z),  Ng  Is  the  order  of  B(z),  M is  the  number 
of  closed  loop  characteristic  roots  out  of  the  unit  circle  and  y is  the 
number  of  open  loop  characteristic  roots  out  of  the  unit  circle.  In 
detail,  the  total  phase  change  of  P(z)  along  r can  be  evaluated  by  the 
summation  of  Individual  phase  changes  along  4(1),  i ■ + I,  ♦ 2,  ...,+^n 
and  semicircles  c(l),  1 » 0,  + 1,  ♦ 2,  . . . , ♦ (n-1),  n.  In  other  words 


Assume  that 


- 5 4et(.)  * av.>  - ,X, 


P(ej“T)  - U(u>)  + jV(co) 


where  U(<o)  and  V(u)  are  real  functions  of  u>.  From  property  (3),  we  have 


N N. 

1(2)  z-P>ane 

'Veto 

\i(i) 

Vc(°) 

JO — H 

c(nM 

l(-n)V 

/(-I) 

yC(-D 

t(-2) 

FIGURE  5-12.  THE  NYQUIST  CONTOUR  FOR  DISCRETE-TIME  CONTROL  SYSTEMS 
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P*(eJwT)  - U(u»)  - jV(u) 

- P(e"Ju,T) 

- U(-u)  + jV(-u) 

Therefore,  following  the  definition  of  the  phase  angle  change,  we 

aVd  " AV-i) 

Substituting  Eq.  (66)  into  Eq.  (63),  yieids 

n n 

A0r  “ 2 E Ae»m  +2^0  n\ 

F 1-1  4(1)  i-Tn-i)  c(,) 

Next,  consider  the  phase  angie  change  along  the  semicircles  c(l). 
that  there  are  Bj  poles  of  P(z)  at  z - e^“t , then 

a0cd)  ■ *V 

Thus 

"•«'»  ‘ " ,-t-o  e' 

n 

where  B ■ E B.  «nd  Is  the  total  number  of  poles  of  P(z)  on 
I— Tn-1)  1 

unit  circle.  Combining  Eqs.  (62),  (67)  and  (69)  yields 

S A6l(l)  " T 12(y-M‘K)  + 
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(65) 

have 

(66) 


(67) 

Assume 


(68) 


(69) 

the 


(70) 


a 


r 


K - nb  ■ na 


Now,  we  can  state  the  generalized  Nyqulst  criterion  for  discrete-time 
control  systems,  as  follows: 


If  the  return  difference  P(z)  of  the  discrete-time  system  has 
the  properties  (1)  through  (3),  the  system  Is  stable  if  and 
only  If  the  total  phase  change  along  the  open  arcs  t(l), 

I ■ 1,  2,  ....  n.  Is  equal  to  y [2(y-K)  + 6]. 

Substituting  H - 0 Into  Eq.  (70),  we  Immediately  prove  the  criterion. 

Program  No.  9 (Appendix  C)  was  written  for  use  on  the  Fourier  Analyzer  to 

calculate  and  plot  the  phase  angle  of  the  return  difference  P(z)  for 
u 

o < u < y~  . 

Example  6.  Consider  the  system  with  transfer  functions 


-zA  - 2z3  + 0.7z2  + O.iz  ♦ A 


(zZ  - l)(zZ  + A) 


H(z)  - 1 


The  return  difference  Is  given  by 


P(z)  - 


z3  + 0.7z2  ♦ 0.1: 


+ 3z  - A 


From  Eqs.  (72)  and  (7*0,  we  have 


6-2,  y - 2,  K - 1 


Therefore,  the  stability  criterion  Is 
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A6  - j [2 (y-K)  + b!  - 2n 


The  phase  plot  of  Eq.  (J1*)  Is  shown  In  Figure  5” 1 3 - The  phase  change  Is 
2ir;  therefore,  the  system  Is  stable. 

Example  7-  Consider  the  discrete-time  system  with  transfer  functions 


G(z)  - 


(z  - 3) (z  + 3)(z  + 1) 


H (z)  - 1 


The  return  difference  Is  given  by 


Examining  Eqs.  (75)  and  (77),  yields 


6-2,  y - 2,  K-0 


Thus,  the  stability  criterion  Is 


A6  - y (2 (y-K)  + B] 


PHASE  PLOT  OF  THE  RETURN  DIFFERENCE  OF  EXAMPLE  6 


The  phase  plot  of  Eq.  (77)  Is  shown  in  Figure  5-14.  The  total  phase 
change  Is 


A0  - A0,  + A02 

- 270°  + 270° 

- 3* 

Therefore,  the  system  is  stable. 

Example  8.  Consider  the  stability  of  the  missile  launching  system  given 
in  Section  5.**.  The  return  difference  is  given  by 

P(z)  - 1 + Dc(z)  G(z) 

Figure  5~ 1 5 shows  the  phase  plot  of  the  return  difference;  the  phase  change 
Is  ~ . From  Eqs.  (45)  and  (49),  we  have 

t 

0*1,  Y-0,  K - 0 

According  to  the  generalized  Nyquist  criterion 


AO  • j l2(y-K)  + B]  - j 


Therefore,  the  system  is  stable. 


Example  9.  Consider  the  discrete-time  system  with  transfer  function 

(78) 

H(z)  - 1 (79) 


G(z)  * -(zTlTlz2-  0.5) 


The  return  difference  Is  given  by 
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P(z) 


(80) 


z + 0.5z  - 0.5 

Examining  Eqs.  (78)  and  (80),  yields 

6-1,  Y-0,  K-0 

Thus  the  stability  criterion  Is  A0  » j . The  phase  plot  Is  shown  In 
Figure  5“ 1 6 ; the  phase  angle  change  Is  - y . Therefore,  the  system  Is 
unstable. 

5.6  REMARKS 

A unified  approach  to  discrete-time  control  system  analysis  via  the 
FFT  algorithm  has  been  established  In  this  chapter.  Even  though  only 
slngle-lnput/slngle-output  systems  were  studied,  It  Is  evident  that  the 
approach  can  be  easily  extended  to  the  multiple  input/output  case. 
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PHASE  PLOT  OF  THE  RETURN  DIFFERENCE  OF  EXAMPLE  8 
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CHAPTER  VI 

DESIGN  OF  DISCRETE-TIME  CONTROL  SYSTEMS 


6.0  INTRODUCTION 

A design  procedure  for  discrete-time  feedback  control  systems  is 
presented  In  this  chapter:  (1)  to  study  the  closed  loop  specification 

by  converting  various  requirements  into  a unified  language  which  Is 
manageable  via  the  minicomputer;  (2)  to  construct  a reference  model,  for 
example,  a required  closed  loop  frequency  response;  (3)  to  search  for 
the  suitable  frequency  response  of  the  compensator  via  the  FFT  algorithm; 
(4)  to  design  a nonrecurslve  filter  transfer  function  for  the  compensator 
using  the  window  method;  and  (5)  to  verify  the  performance  of  the  com- 
pensated system.  Several  design  examples  including  type  "0"  and  type 
”1"  systems  are  developed  for  illustration. 

6.1  FREQUENCY  DOMAIN  APPROACH  FOR  DISCRETE-TIME  CONTROL  SYSTEM  DESIGN 

Consider  a basic  configuration  of  a direct  digital  control  (DDC)  sys- 
tem as  shown  in  Figure  6-1.  Gp(s)  is  the  plant;  H(s)  Is  the  feedback 
device;  and  Gq(s)  Is  the  holding  circuit  or  D/A  converter.  Usually  the 
zero  order  hold  is  taken  to  be 

1 - e"5T 

«„<*>  • •’-r—  "> 

where  T is  the  sampling  period.  0 (z)  Is  the  digital  compensator  and  Z 

, c 

is  the  sampler  or  A/D  converter. 

Using  the  block  diagram  operations  described  In  Chapter  V,  a simpli- 
fied block  diagram  In  the  z-domaln  is  obtained  and  Is  shown  In  Figure  6-2, 
where 


G,(z)  - Z[Gp(s)  Gq(s)] 


G2(z)  - Z[Gq(s)  Gp(s)  H(s)] 


Therefore,  the  closed  loop  transfer  function  Is  given  by 


rM  c(z)  °e(*>  M‘> 

F(z)  ’nr  r*  oc(z1c2T; 


The  design  problem  Is  to  find  the  compensator  D £(z)  such  that  the  charac- 
teristics of  the  compensated  system  match  the  required  specifications. 

If  Gj(z)  and  G2(z)  are  given,  from  Eq.  (b)  we  have 


Dc(z)  - 


In  other  words,  If  the  closed  loop  transfer  function  F(z)  Is  specified, 

the  compensator  D (z)  Is  given  by  Eq.  (5).  However,  there  are  some  prob- 
c 

lem  areas: 


(1)  The  analytic  form  of  the  specific  closed  loop  transfer 
function  is  necessary. 


(2)  The  compensator  D (z)  determined  by  Eq.  (5)  Is  In 

c 

recursive  form  and  since  the  characteristics  of  a 
recursive  filter  are  very  sensitive  to  Its  coefficients, 
high  accuracy  In  determining  coefficients  Is  necessary. 
Therefore  the  microcomputer,  which  usually  has  a data 
word  length  of  less  than  16  bits,  is  not  always  suitable 
to  be  employed  as  a recursive  digital  compensator. 

(3)  The  compensator  determined  by  Eq.  (5)  may  be  unstable. 


In  order  to  avoid  the  disadvantages  of  the  design  method  in  z-domain, 
we  will  develop  a frequency  domain  approach.  The  frequency  response  of 
Dc(z)  is  given  by 


,(e>T)  - 


i,(cj“T)  - 


F(eJtoT) 
G2(eJuT 


) F (eJuT) 


(6) 


If  Gj  (z)  and  G2(z)  are  given,  the  frequency  response  G^  (e^“T)  and  G2(e-*wT) 
can  be  found  using  the  FFT  algorithm  as  described  in  Chapter  V.  The 
required  information  concerning  the  closed  loop  characteristics  is  the 
frequency  response  F(e"*w^).  Therefore,  F(z)  can  be  specified  either  in 
analytic  form  or  by  the  time  or  frequency  domain  characteristics  from 
which  an  approximate  frequency  response  can  be  found.  If  we  take  the 
Inverse  Fourier  Transform  of  D (e^“^) , we  can  easily  find  the  non- 
recursive filter  form  for  D (z).  In  other  words,  the  window  method  for 
nonrecurslve  digital  filter  design  as  given  in  Chapter  IV  can  be  used 
to  design  a discrete-time  control  system. 


6.2  SPECIFYING  THE  FREQUENCY  RESPONSE  OF  CLOSED  LOOP  SYSTEMS 


Many  approaches  are  available  for  specifying  the  closed  loop  fre- 
quency response  to  match  the  required  specifications.  Two  of  these 
approaches  are  discussed  in  the  ensuing  paragraphs. 


(I)  The  Required  Specification  is  Given  In  The  Continuous  Time 
Domain  or  The  Frequency  Domain 

When  the  output  of  the  discrete-time  control  system  Is  con- 
tinuous In  time,  we  can  specify  the  closed  loop  transfer 
function  In  the  s-domatn  according  to  the  required  specifi- 
cations and  then  find  its  transfer  function  in  z-domain. 

For  example.  If  we  need  the  closed  loop  system  which  has 
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C - 0.7 


“d  “ 5 
Ky  - 20 


then,  by  the  use  of  the  multi -dimensional  Newton-Raphson 

[273 

method  the  required  closed  loop  transfer  function  can 
be  determined  and  is  given  by 


F (s)  - M50s  ♦ 12,674 
1 s2  + 4.984s  + 12.674 


(7) 


Assuming  the  sampling  period  of  the  system  is  T *■  0.1  seconds; 
the  z-transform  of  F(s)  Is  given  by 


F(z)  - 4.3502  - 3.l40z 


1.5087z  + 0.6075 


(8) 


Using  the  FFT  algorithm,  the  frequency  response  of  the  closed 
loop  transfer  function  can  be  easily  found. 

(2)  The  Specification  is  Given  by  The  Discrete  Transient  Response 

If  the  impulse  response  of  the  closed  loop  is  specified,  the 
frequency  response  can  be  obtained  by  taking  the  inverse 
Fourier  Transform  of  the  tmpulse  response.  Sometimes,  the 
step  response,  instead  of  the  impulse  response,  of  the  closed 
loop  Is  specified.  In  this  case,  we  need  to  find  the  impulse 
response  first.  Since  the  step  response  is 


Cu(z)  - F(z) 

the  impulse  response  is  given  by 


(9) 
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C4w  ■ cu(z) 


2 - 1 


z 

•cu(z)0-z"1)  00) 

In  the  time  domain,  we  have 

c4(t)  - Cu(t)  - Cu(t-T)  (II) 

The  Fourier  Analyzer  contains  a differentiator  which  is  actually 
an  operator  with  transfer  function  1 - z *;  therefore,  the 
impulse  response  can  be  found  from  the  step  response  using  the 
differentiator.  Once  the  impulse  response  is  obtained,  the 
frequency  response  of  the  closed  loop  can  be  found  easily. 
Program  No.  10  (Appendix  C)  can  be  used  to  find  the  frequency 
response  from  the  step  response. 


For  example,  the  required  unit  step  response  is  listed  in  Table  6-1 
and  the  sampling  period  Is  T • 0.1  seconds.  The  transient  response  for 
0 _<  t <^5.0  is  plotted  in  Figure  6-3(a).  Using  Program  No.  10  (Appendix 
C),  the  impulse  response  was  obtained  and  is  shown  in  Figure  6-3(b)  and 
listed  In  Table  6-1.  The  frequency  response  Is  shown  in  Figures  6-3 (c) 
and  6-3 (d).  — _ 


TABLE  6-1 

REQUIRED  UNIT  STEP  AND  IMPULSE  RESPONSES 


t 

cu(t) 

c6(t) 

0 

0 

0.1 

0.280 

0.2 

0.570 

0.3 

0.200 

o.<» 

0.050 

0.5 

1.070 

-0.030 

0.6 

1.030 

-0.01*0 

0.7 

0.990 

-0.01*0 

0.8 

0.980 

-0.010 

t 

Cu(t) 

c6(t) 

0.985 

0.005 

mm 

0.995 

0.010 

i.i 

0.015 

1.2 

1.005 

-0.005 

1.3 

0.998 

-0.007 

1.1* 

0.999 

0.001 

1.5 

1.000 

0 

a 

1.000 

0 

a 

00 

1.000 

0 
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6.3  THE  STEADY  STATE  ERROR  AT  THE  SAMPLING  INSTANTS 

Consider  the  unit  feedback  system  as  shown  in  Figure  6-4.  The  sys 
tern  error  is  defined  as 


E(z)  - R(z)  - C(z) 


R(z) 

ttttff 


(12) 


where  R(z)  Is  the  input  and  C * «)  Is  the  output.  From  the  Final  Value 
theorem,  the  steady  state  error  Is  given  by 


*ss  " ' im 


- Urn  E(z) 

z-1  Z 


» I Im 

z-1 


1)  E(z) 

rcdTT 


(13) 


FIGURE  6-4.  THE  UNIT  FEEDBACK  SYSTEM 

The  following  three  kinds  of  steady  state  errors  are  most  Important  In 
considering  the  closed  loop  specifications: 


14! 


(1)  Position  errors. 


The  position  error  e^  is  defined  as  the  steady  state  error 
when  the  input  is  the  unit  step  function 


From  Eq.  (13).  we  have 

®P  " Ji"  ttgtit  " r ♦1g(i>  " IT 

where 


O'*) 


(15) 


Kp  - 1 + G ( 1 ) (16) 

is  called  the  position  error  constant.  For  type  "0"  systems, 
the  position  error  ep  is  finite;  for  type  "1"  or  higher  sys- 
tems, the  position  error  e Is  zero. 

P 

(2)  Velocity  errors. 

The  velocity  error  ev  Is  defined  as  the  steady  state  error 
when  the  system  Is  excited  by  the  unit  ramp  signal 

M*>  (17) 

(z  - I)2 

Thus,  from  Eq.  (13)  Me  have 

•v  " 77-  mf  vd(z)j 

1 


(18) 

(19) 


1<*2 


where 


K - * 1 lm  (z  - l)  G(z)  (20) 

V 1 z-1 

Is  called  the  velocity  error  constant.  For  type  "0"  systems 
ev  ■ *5  for  type  "1"  systems  ey  Is  finite;  and  for  type  "2"  or 
higher  systems  ey  ■ 0. 


(3)  Acceleration  errors. 

The  acceleration  error  e#  Is  defined  as  the  steady  state 
error  when  the  system  Is  excited  by  the  unit  parabolic  sig- 
nal 

r(t)-It2  (21) 


or 


A(z) 


, i2(*  ♦ »>* 

2 (z  - l)2 


Therefore,  from  Eq.  (13)  we  have 

T2 


e„  ■ 1 Im r 

a z-»l  (z  - 1)2ll  + G(z)l 


where 


K ■ ir  11m  (z  - l)2  G(z) 

a r z-i 


(22) 


(23) 


(24) 
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Is  called  the  acceleration  error  constant.  For  type  "0"  and 

type  "1"  systems,  e ■ «•;  for  type  "2"  systems  e Is  finite; 

o a 

for  type  "3"  or  higher  systems  ea  * 0. 


If  the  open  loop  transfer  function  G(z)  Is  given.  It  Is  straight- 
forward to  find  the  error  constants  K , K and  K . If  only  the  closed 

p v a 

loop  transfer  function  F(z)  Is  given,  we  have  to  find  the  open  loop 
transfer  function 


G(z) 


F(z) 

't~fTzT 


(25) 


and  then  find  the  error  constants.  However,  If  only  the  frequency  response 
of  the  closed  loop  Is  specified  as  described  In  Section  6.2,  we  do  not 
have  the  analytic  form  for  F(z),  and  must  find  the  error  constants  from 
the  frequency  response. 


Substituting  Eq.  (25)  Into  Eq.  (16),  yields 


k ■ 1 im 
p z-1 


[’ + r^rta] 


(26) 


Since  z ■ e^u^,  z-*l  corresponds  to  u-> 0.  Thus,  In  the  frequency  domain 
Eq.  (26)  becomes 


1 


K ■ llm  i-  y 

P ur*0  1 - F(eJuT) 


(27) 


Similarly,  from  Eqs.  (20)  and  (2*0 , we  have 


k . ii. 

V wr+0  T [l  “ F(eJwT)] 


(28) 


and 
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(Mma 


K HeJ“T)(«J“T  -l.) 

8 o)-*0  T2  [l  - F(eJuT)] 


(29) 


Most  servo  systems  or  autopilot  systems  are  type  "1",  and  the  velocity 
error  is  important;  therefore,  we  will  discuss  the  velocity  error  con- 
stant in  detail.  For  the  type  1 system,  the  open  loop  transfer  func- 
tion can  be  rewritten  as 


G(z)  - 


X(z) 


(30) 


where  X(l)  is  nonzero  finite.  Thus,  the  closed  loop  transfer  function 
Is  given  by 


-/  \ G(z)  _ X (z) 
F(z)  * 1 ' + G( TT  (z  - )) 


(31) 


and 


F(l)  - I 


(32) 


In  other  words,  F(z)  - 1 has  a factor  (z  - 1).  From  Eq.  (31),  we  have 


dFW  u-o^-xw 


(33) 


Therefore,  we  have 


F'(,)  " ' xTTT  (3 

Since  X (1 ) is  nonzero  finite,  we  conclude  that  F'(!)  - 0 and  z - 1 Is 
a simple  root  of  F(z)  - 1 ■ 0 or 


1*5 


Is  nonzero  finite  as  z-*l . 


Y(*)  £ ^4-y-1 


(35) 


If  we  define 


D(z) 


F(z) 


(36) 


then  z ■ 1 is  the  simple  root  of  F(z)  -1*0,  and  the  denominator  of 
Eq.  (36)  Is  nonzero  finite  as  z+1 . From  Eq.  (28)  we  have 


Kv  - I 


(37) 


Another  approach  to  finding  Kv  Is  to  first  find  the  error  response  E(z). 
From  Eq.  (12),  the  error  response  of  the  unit  ramp  Input  is 


E(z) 


Tz 


(z 


- i r 

TT7J— 

I - F(zV 


1 + 


(38) 


Since  z~ F(z)  Is  the  unit  step  response  of  the  closed  loop,  y - 
z ~z-  | F(z)  is  the  error  response  of  the  unit  step  Input.  If  the  closed 
loop  characteristic  Is  specified  by  Its  unit  step  response,  the  error 
response  of  the  unit  step  input  in  the  time  domain 


1l*6 


can  be  found  easily.  Then,  the  error  response  of  the  unit  ramp  in  the 
time  domain 


er(t)  - z”  [E(z)] 

can  be  obtained  by  integrating  ( ■*  ^ ) and  the  right  shifting  (Z  ') 
of  Eu ( t ) . Once  e r ( t ) is  found,  the  velocity  error  Is 


(AO) 


e - 11m  [t  - er(t)] 
t-**° 


(Al ) 


and 


K -JL 

v e 

v 


(A2) 


The  entire  procedure  for  finding  the  frequency  response  and  analyzing 
the  velocity  error  from  the  specified  step  response  of  the  closed  loop 
is  included  in  Program  No.  11  (Appendix  C) . 

I 

For  the  system  with  closed  loop  step  response  shown  In  Figure  6-3(a), 
the  frequency  response  of  D(z)  [Eq.  (36)1  Is  shown  in  Figure  6-A(a). 

From  Figure  6*A(a),  we  have  D(e^uT)  ■ 0.60A  when  u - 0;  therefore. 


DO)  _ 


K - 6.0A 

v T 


(A3) 


The  error  response  of  the  unit  ramp  Input  is  shown  In  Figure  6-A(b). 
The  steady  state  error  Is  approximately  0.165A;  thus, 


K - i-  - 6.0A595 
v ev 


and  the  results  of  Eqs.  (A3)  and  (AA)  are  consistent. 


(AA) 


1A7 


J 


6.1*  CALCULATION  OF  THE  FREQUENCY  RESPONSE  OF  REQUIRED  COMPENSATOR 

The  basic  formula  for  finding  the  frequency  response  of  the  required 
compensator  is  given  by  Eq.  (6).  For  computing  advantages,  we  need  to 
study  the  following  system  categories: 

(1)  Type  "0"  systems. 


For  the  type  "0"  system,  no  pole  occurs  at  z * 1 for  both 
G| (z)  and  G2 (z) ; therefore,  Eq.  (6)  can  be  used  without 
difficulty.  In  the  unit  feedback  system  as  shown  In 
Figure  6-1*,  G(z)  - Gj(z)  ■ G2(z);  therefore,  we  have 


D (eJ-T,  . F(eJtoT) 

c G(ej“T)[l  - F(eJuT)3 


When  u ■ 0,  or  e-*uT  ■ 1 , we  have 

Dc(,)  “ emi^-Vror 
From  Eq.  (27)  we  have 


p " nm 


and 


F ( 1 ) -1+J- 

P 


Thus,  we  have 


1 + K 

°c(,)  "-gITT* 


(1*5) 


(W) 


(1*7) 


(1*8) 


(1*9) 


1 1*8 


When  a small  position  error  or  a large  position  error 
constant  (K^)  Is  required,  Eq.  (1*9)  indicates  that  the 
dc  gain  of  the  frequency  response  of  the  compensator 
should  be  as  large  as  possible. 

(2)  Type  "1"  systems. 

For  type  "1"  systems,  G^z)  and  G2(z)  become  infinite 
when  u-*0  or  z-*-l . As  shown  in  the  last  sect  ion,  [1  - F(z) ] 
has  a factor  (z  - 1).  Thus,  in  the  unit  feedback  system, 
G(z)  ■ Gj(z)  = G2(z);  Eq.  (1*5)  gives  an  undefined  value 
for  D£(z)  when  w-*0.  In  order  to  avoid  this  ambiguity, 
we  must  modify  the  formula  for  the  frequency  response 
of  the  compensator. 

Assume  that 


G(z>  -Htc„w 


(50) 


when  Gp (z)  has  no  pole  or  zero  at  z = 1.  From  Eq.  (1*5) , 
we  have 


Substituting  Eq.  (36)  Into  Eq.  (51),  yields 


»C<*> 


0(2) 


(52) 


Since  the  frequency  response  of  D(z)  is  well-defined  for 
0 <_u>  <_ws,  the  frequency  response  of  the  required  com- 
pensator Is  given  by 
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Dc(eJ“T) 


0(eJ“T) 

Gn(«J“T) 


From  Eq.  (37),  we  see  that 


(53) 


0(1)  - TKy  (5*> 

When  u * 0,  we  have 

TK 

°co)  • infr  (55) 

n 

Eq.  (55)  Indicates  that  the  dc  gain  of  the  frequency  response 

of  the  compensator  Is  proportional  to  the  velocity  error 

constant  K . 

v 

6.5  GENERAL  PROCEDURE  FOR  DESIGNING  NONRECURSIVE  COMPENSATORS  BY  THE 
WINDOW  METHOD 


Once  the  frequency  response  of  the  required  compensator  is  found, 
the  frequency  response  can  be  expanded  as  a Fourier  series  In  the 
frequency  domain  and  then  truncated  through  the  use  of  windows  to  con- 
struct a nonrecursive  compensator  of  the  proper  order.  The  basic 
principle  is  the  same  as  the  window  method  for  nonrecursive  digital 
filter  design. 

As  pointed  out  In  Chapters  II  and  III,  the  effects  of  windowing 
are:  (1)  to  reduce  or  to  eliminate  the  Gibbs'  phenomenon  in  the  fre- 

quency domain  and  (2)  to  improve  the  damping  rate  of  the  side  lobes  of 
the  rectangular  window  which  corresponds  to  the  trucatlon  of  the  Fourier 
series.  In  most  cases,  we  do  not  have  a point  of  discontinuity  in  the 
frequency  response  of  the  compensator;  therefore,  the  principal  effect 
of  windowing  in  nonrecursive  compensator  design  is  the  latter.  In 
fact,  this  effect  smooths  the  frequency  response  of  the  compensator. 
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The  general  procedure  for  nonrecursive  compensator  design  by  the 
window  method  can  be  outlined  as  follows: 

(1)  Specify  the  closed  loop  characteristics;  such  as  the 
overshoot,  the  rise  time,  the  settling  time,  the 
velocity  error,  the  position  error,  etc. 

(2)  Find  the  frequency  response  of  the  closed  loop  matching 
the  required  specifications. 

(3)  Calculate  the  frequency  response  of  the  required  compen- 
sator from  the  frequency  responses  of  the  open  loop  and 
the  closed  loop. 

(*•)  Construct  the  transfer  function  of  the  compensator  from 
Its  frequency  response  by  the  window  method  and  deter- 
mine the  compensator  in  nonrecurslve  form. 

(5)  Check  the  closed  loop  characteristics  of  the  compensated 
system;  such  as  the  frequency  response,  the  impulse 
response,  the  step  response,  etc. 

The  whole  procedure  can  be  performed  using  the  FFT  algorithm.  Following 
this  procedure,  the  Fourier  Analyzer  becomes  a power  tool  for  discrete- 
time  control  system  design. 

Since  there  are  some  differences  In  computation  of  the  frequency 
responses  of  the  compensator  for  type  "0"  and  type  “1"  systems,  two 
separate  programs  were  prepared  for  the  design  of  each  type.  Program 
No's  12  and  13  (Appendix  C)  were  written  for  type  "0"  and  type  "1" 
system  designs  respectively. 
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(c)  FREQUENCY  RESPONSE  OF 
THE  REQUIRED  COMPENSATOR 


(b)  REQUIREO  CLOSED  LOOP 
FREQUENCY  RESPONSE 


FIGURE  6-7.  DESIGN  MODELS  OF  EXAMPLE  1 


Djtz)  which  Is  obtained  using  the  Hanning  window,  the  corresponding  sys- 
tem frequency  response,  and  the  unit  step  response  of  the  compensated 
system  with  compensator  D2(z),  respectively.  From  Figure  6-9,  the  rise 
time  Is  approximately  3T;  the  overshoot  is  approximately  10%,  and  the 
position  error  constant  Is  approximately  A. A.  Figures  6- 10(a) , (b)  and 
(c)  show  the  frequency  response  of  Dj(z)  which  is  obtained  using  the 
trapezoidal  window,  the  corresponding  system  frequency  response,  and  the 
unit  step  response  of  the  compensated  system,  respectively.  From 
Figure  6-10,  the  rise  time  Is  approximately  3T;  the  overshoot  is  approxi- 
mately 12%;  and  the  position  error  constant  is  approximately  A. 9. 


From  this  example,  we  see  that  the  effect  of  windowing  is  remarkable. 
Also,  we  note  that  the  trapezoidal  window  with  parameter  b = 0.A  performs 
as  well  as  the  Hanning  window. 

Example  2.  Consider  the  system  as  shown  In  Figure  6-11.  The 
required  specifications  are  as  follows: 


(1)  The  rise  time  Is  less  than  3T,  T ■ 1 second, 

(2)  The  overshoot  of  the  unit  step  response  is  less  than  30%,  and 

(3)  The  position  error  constant  Kp  « 10. 

The  open  loop  transfer  function  In  the  z-domaln  Is  given  by 


G(z)  * z fs~  httvtt 


0.2326z 

- 0.5032z  + 0.0A978 


(57) 


The  frequency  response  of  G(z)  is  shown  in  Figure  6-12. 

Using  Figure  6-!3(a)  as  a model  of  the  unit  step  response  for  the 
compensated  closed  loop  system,  the  required  frequency  response  is  shown 
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(a)  FREQUENCY  RESPONSE  OF  (b)  FREQUENCY  RESPONSE  OF 

THE  COMPENSATOR  D2(z)  THE  COMPENSATED  SYSTEM 


(c)  UNIT  STEP  RESPONSE  OF  THE  COMPENSATED  SYSTEM 
FIGURE  6-9.  PERFORMANCE  OF  THE  S”S* EM  COMPENSATED  WITH  D,(z) 


In  Figure  6- 1 3 (b) - The  frequency  response  of  the  compensator  Is  shown 
In  Figure  6- 1 3 (c) . The  coefficients  of  nonrecurslve  compensators  obtained 
using  various  windows  are  listed  In  Table  6-3. 


Figures  6-H(a),  (b)  and  (c)  show  the  frequency  response  of  Dj  (z) 
obtained  by  using  the  rectangular  window,  the  corresponding  system  fre- 
quency response,  and  the  unit  step  response  of  the  compensated  closed 
loop  system,  respectively.  By  observation,  too  many  ripples  appear  in 
the  closed  loop  frequency  response  and  too  many  wavelets  in  the  transient 
response;  therefore,  this  compensator  is  not  acceptable. 

Figures  6-l5(a),  (b)  and  (c)  show  the  frequency  response  of  DjCz) 
obtained  using  the  Hanning  window,  the  corresponding  system  frequency 
response  and  the  unit  step  response  of  the  compensated  system  respec- 
tively. From  Figure  6-15,  the  rise  time  is  approximately  2T;  the  over- 
shoot Is  approximately  28%;  and  the  position  error  constant  Kp  ® 9. 

Figures  6-i6(a),  (b)  and  (c)  show  the  frequency  response  of  D^(z) 
obtained  using  the  trapezoidal  window  with  parameter  b « O.b,  the 
corresponding  system  frequency  response,  and  the  unit  step  response  of 
the  comepnsated  system,  respectively.  Figure  Figure  6-16,  the  rise  time 
is  approximately  2T;  the  overshoot  is  approximately  30%;  and  the  position 

error  constant  K 9 10. 

P 

The  results  of  these  two  examples  were  determined  using  Program  No. 
12  (Appendix  C). 

6.7  TYPE  "1"  SYSTEM  DESIGN 

[22] 

Example  3.  Consider  the  infrared  tracking  system  as  shown  In 
Figure  6-17.  The  Input  r(t)  Is  the  target  1 Ine-of-slght;  the  output 
c(t)  is  the  tracker  1 Ine-of-slght;  D£(z)  Is  the  digital  compensator  and 
Gq(s)  Is  the  zero  order  hold  given  by 
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(c)  FREQUENCY  RESPONSE  OF 
THE  REQUIRED  COMPENSATOR 


FIGURE  6-13.  DESIGN  MODELS  OF  EXAMPLE  2 


,tr  :•  i 


I 


TABLE  6-3 

COMPENSATORS  FOR  EXAMPLE  2 


LEGEND 

D,  (z)  - ^ d,  2’-*,  for  1-1,2,  3. 

1 J-0  J 

Dj(z)  - Designed  by  using  the  rectangular  window. 

D^tz)  - Designed  by  using  the  Hanning  window. 

D„(z)  - Designed  by  using  the  trapezoidal  window,  b 
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(a)  FREQUENCY  RESPONSE  OF 
THE  COMPENSATOR  D? (z) 


(b)  FREQUENCY  RESPONSE  OF 
THE  COMPENSATED  SYSTEM 


! 


(c)  UNIT  STEP  RESPONSE  OF  THE  COMPENSATED  SYSTEM 


FIGURE  6-1  A.  PERFORMANCE  OF  THE  SYSTEM  COMPENSATED  WITH  D, (z) 
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(a)  FREQUENCY  RESPONSE  OF 
THE  COMPENSATOR  D2(z) 


(b)  FREQUENCY  RESPONSE  OF 
THE  COMPENSATED  SYSTEM 


(c)  UN^T  STEP  RESPONSE  OF  THE  COMPENSATED  SYSTEM 


FIGURE  6-15.  PERFORMANCE  OF  THE  SYSTEM  COMPENSATED  WITH  D2(z) 
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(a)  TRACKING  SYSTEM 


(b)  SERVOMOTOR  LOOP 

FIGURE  6-17.  INFRARED  TRACKING  SYSTEM 


168 


1 "ST 

60(s)  - — (58) 

The  transfer  function  of  the  servomotor  loop  Is 

GP(S)  " STs  V 3.5)'  (59) 

The  required  specifications  of  the  IR  tracking  system  are  as  follows: 

(1)  The  overshoot  is  less  than  25%, 

(2)  The  rise  time  is  less  than  3T, 

(3)  The  settling  time  is  less  than  15T, 

(i»)  The  velocity  error  constant  is  greater  than  10  sec  and 
(5)  The  bandwidth  is  greater  than  6.28  rad/sec. 

t 

According  to  Eq.  (5),  we  must  choose  the  sampling  frequency,  u>s  > 2wb  = **it. 
Thus,  the  sampling  frequency  should  be  greater  than  2 Hz.  Let  us  choose 
fs  - 10  Hz,  or  a sampling  period  of  T ■ 0.1  seconds. 


The  open  loop  transfer  function.  In  the  z-domain,  of  the  tracking 
loop  Is  given  by 


G(z)  - Z[G0(s)  Gp(s) ] 

0.00i»A6(z  + 0.8903) 
“ (z  - IMz  - 0.70^71 


(60) 


Using  Figure  6- 1 8 (a)  as  the  model  of  the  unit  step  response  for  the 
compensated  closed  loop  matching  specifications  (1)  through  (I*),  the 
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corresponding  closed  loop  frequency  response  is  shown  In  Figure  6- 18(b) . 
The  frequency  response  of  the  required  compensator  Is  shown  In  Figure 
6- 1 8 (c) . The  coefficients  of  nonrecursive  compensators  obtained  using 
various  windows  are  listed  in  Table  6-4. 

Figures  6- 19(a) , (b)  and  (c)  show  the  frequency  response  of  (z) 
obtained  using  the  rectangular  window,  the  corresponding  system  frequency 
response  and  the  unit  step  response  of  the  compensated  system  respec- 
tively. By  observation,  too  many  ripples  appear  in  the  closed  loop 
frequency  response  and  too  many  small  wavelets  in  the  transient  response; 
therefore,  this  compensator  is  not  suitable. 

Figures  6-20 (a),  (b)  and  (c)  show  the  frequency  response  of  D2 (z) 
obtained  using  the  Hanning  window,  the  corresponding  system  frequency 
response  and  the  unit  step  response  of  the  compensated  system  respec- 
tively. From  Figure  6-20  we  find:  (1)  the  overshoot  Is  20$,  (2)  the 

rise  time  Is  2T,  (3)  the  settling  time  is  10T,  and  (4)  the  velocity 
error  constant  is  - 9.97. 

t 

Figures  6-21  (a),  (b)  and  (c)  show  the  frequency  response  of  D^(z) 
obtained  using  the  trapezoidal  window  with  the  parameter  b = 0.4,  the 
corresponding  system  frequency  response  and  the  unit  step  response  of 
the  compensated  system  respectively.  From  Figure  6-21  we  find:  (1)  the 

overshoot  Is  22$,  (2)  the  rise  time  is  2T,  (3)  the  settling  time  is  10T, 
and  (4)  the  velocity  error  constant  Is  ■ 11.2. 

Comparing  Figure  6-21  with  Figure  6-20,  we  note  that  the  trapezoidal 
window  performs  even  better  than  the  Hanning  window  and  yields  a larger 
velocity  error  constant. 
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(b)  REQUIRED  CLOSED  LOOP  (c)  FREQUENCY  RESPONSE  OF 

FREQUENCY  RESPONSE  THE  REQUIRED  COMPENSATOR 


FIGURE  6-18.  DESIGN  MODELS  OF  EXAMPLE  3 


TABLE  6-4 


COMPENSATORS  FOR  EXAMPLE  3 


COEFFICIENTS 

Dj  (z) 

0,  (z) 

Oj(z) 

d0 

66.8 

66.8 

66.8 

dl 

-32.8 

-32.4 

-32.8 

d2 

-40.7 

-39.2 

-40.7 

"3  " 

76.7 

70.1 

76.7 

d4 

-80.4 

-68.6 

-74.2 

d5 

79.2 

61.4 

67.0 

d6 

-67.9 

-46.9 

-52.0 

d7 

63.0 

37.5 

43.6 

d8 

-54.2 

-27.1 

-33.4 

d9 

52.  1 

20.9 

28.0 

dl0 

-43.9 

-13.5 

-20.2 

dll 

40.8 

9.0 

15.7 

d12 

-32.6 

- 4.8 

-10.0 

dl3 

30.0 

2.5 

6.9 

dl4 

-22.85 

- 0.9 

- 3.5 

d15 

21.4 

0.2 

1.7 

LEGEND 


°I(z)  “ z'J,  for  1-1,2,  3. 

Dj (z)  - Designed  by  using  the  rectangular  window. 

D^z)  - Designed  by  using  the  Hanning  window. 

D^(z)  - Designed  by  using  the  trapezoidal  window,  b - 0.4. 
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THE  COMPENSATOR  D,  (z)  THE  COMPENSATED  SYSTEM 


(c)  UNIT  STEP  RESPONSE  OF  THE  COMPENSATED  SYSTEM 


FIGURE  6-19.  PERFORMANCE  OF  THE  SYSTEM  COMPENSATED  WITH  D,(z) 
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(a)  FREQUENCY  RESPONSE  OF  (b)  FREQUENCY  RESPONSE  OF 

THE  COMPENSATOR  D2(z)  THE  COMPENSATED  SYSTEM 


(c)  UNIT  STEP  RESPONSE  OF  THE  COMPENSATED  SVSTEM 
FIGURE  6-20.  PERFORMANCE  OF  THE  SYSTEM  COMPENSATED  WITH  D£ (2) 


17^ 


(a)  FREQUENCY  RESPONSE  OF  (b)  FREQUENCY  RESPONSE  OF 

THE  COMPENSATOR  D^z)  THE  COMPENSATED  SYSTEM 


(c)  UNIT  STEP  RESPONSE  OF  THE  COMPENSATED  SYSTEM 
FIGURE  6-21.  PERFORMANCE  OF  THE  SYSTEM  COMPENSATED  WITH  D„(z) 


Example  4.  Consider  the  gun  drive  control  system  as  shown  in 
Figure  6-22,  In  which  the  dynamic  transfer  function  of  gun  drive  is 

Gp(s)  " 7(7  + 3.1)  (s  + 18)  (6,) 

The  required  specifications  are  the  same  as  those  of  Example  3.  Let  us 
choose  a sampling  period  of  T * 0.1  seconds.  The  open  loop  transfer 
function  in  the  z-domain  is 

G(z)  - Z[Gq(s)  Gp(s) ] 

_ 0.005604z(z  - 0. 1A89)  (z  =>  2.342) 

" (z  - 1)  (z  - 0.7047Hz  - 0.1653) 

Using  the  same  model  of  the  unit  step  response  as  shown  in  Figure  6-l8(a), 
the  frequency  response  of  the  required  compensator  Is  shown  in  Figure 
6-23.  The  coefficients  of  nonrecursive  compensators  obtained  using 
various  windows  are  listed  In  Table  6-5. 

Figures  6-24(a),  (b)  and  (c)  show  the  frequency  response  of  Dj(z) 
obtained  using  the  rectangular  window,  the  corresponding  system  frequency 
response  and  the  unit  step  response  of  the  compensated  system  respec- 
tively. From  the  unit  step  response,  we  see  that  the  result  almost  • 
matches  the  required  specifications.  The  velocity  error  constant  is  4 
«v  ■ 18.9.  However,  since  some  ripples  *»re  present  In  the  frecji^^y 
response  of  the  compensated  closed  loop,  the  step  response  has  a small 
wavelet. 

Figures  6-25(a),  (b)  and  (c)  show  the  frequency  response  of  D2(z) 
obtained  using  the  Hanning  window,  the  corresponding  system  frequency 
response  and  the  unit  step  response  of  the  compensated  system  respec- 
tively. From  Figure  6-25,  we  find:  (1)  the  overshoot  is  approximately 

20%,  (2)  the  rise  time  is  2T,  (3)  the  settling  time  Is  7T,  and  (4)  the 
velocity  error  constant  is  14.0.  The  unit  step  response  Is  smoother 
than  that  In  Figure  6-24. 
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TABLE  6-5 

COMPENSATORS  FOR  EXAMPLE  4 


COEFFICIENTS 

D,  (z) 

D,  (z) 

D2(z) 

Dj(z) 

d0 

15.36 

15.36 

15.36 

dl 

17.51 

17.31 

17-51 

d2 

-36.70 

-35.30 

-36.70 

d3 

26.36 

24.08 

26.36 

d4 

-11.64 

- 9.93 

-10.74 

d5 

4.93 

3.82 

4.17 

d6 

- 0.59 

- 0.41 

- 0.46 

d7 

- 0.82 

0.49 

0.57 

d8 

- 0.07 

- 0.04 

- 0.05 

d9 

0.96 

0.39 

0.52 

d10 

0.50 

0.16 

0.23 

dil 

- 0.17 

- 0.04 

- 0.07 

dl2 

1.36 

0.20 

0.42 

d»3 

- 0.65 

- 0.06 

- 0.16 

dl4 

1.47 

0.06 

0.23 

d15 

- 0.56 

- 0.01 

- 0.05 

LEGEND 


0, (*)  - £ d.  z"J,  for  I - 1,  2,  3. 

J-0  J 

Dj (z)  - Designed  by  using  the  rectangular  window. 

D2(z)  - Designed  by  using  the  Hanning  window. 

Dj(z)  - Designed  by  using  the  trapezoidal  window,  b - 0.4. 
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(c)  UNIT  STEP  RESPONSE  OF  THE  COMPENSATED  SYSTEM 


FIGURE  6-21*.  PERFORMANCE  OF  THE  SYSTEM  COMPENSATED  WITH  D, (*) 
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(a)  FREQUENCY  RESPONSE  OF  (b)  FREQUENCY  RESPONSE  OF 

THE  COMPENSATOR  D2 (r)  THE  COMPENSATED  SYSTEM 


(c)  UNIT  STEP  RESPONSE  OF  THE  COMPENSATED  SYSTEM 
FIGURE  6-25.  PERFORMANCE  OF  THE  SYSTEM  COMPENSATED  WITH  D„(z) 


Figures  6-26(a),  (b)  and  (c)  show  the  frequency  response  of  0^(z) 
obtained  using  the  trapezoidal  window  with  the  parameter  b - 0.4,  the 
corresponding  system  frequency  response  and  the  unit  step  response  of 
the  compensated  system  respectively.  From  Figure  6-26,  we  find: 

(I)  the  overshoot  is  approximately  21%,  (2)  the  rise  time  is  2T,  (3)  the 
settling  time  Is  7T,  and  (4)  the  velocity  error  constant  is  15-0.  Also, 
the  unit  step  response  is  smoother  than  that  obtained  by  using  the 
rectangular  window.  Comparing  with  Figure  6-25,  the  step  response  of 
the  compensated  system,  with  the  compensator  obtained  using  a trapezoi- 
dal window  with  b ■ 0.4,  is  approximately  the  same  as  that  found  using 
the  Hanning  window. 

6.8  REMARKS 

The  compensators  designed  in  this  chapter  are  associated  with  direct 
digital  controllers  which  can  be  easily  realized  using  on-line  minicomputer 
or  microcomputer.  Because  of  the  availability  of  the  Fast  Fourier  Trans- 
form algorithm,  the  design  procedure  Is  easier,  simpler,  and  more 
accurate  than  other  existing  methods. 
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(e)  UNIT  STEP  RESPONSE  OF  THE  COMPENSATED  SYSTEM 
FIGURE  6-26.  PERFORMANCE  OF  THE  SYSTEM  COMPENSATED  WITH  D,(x) 

182 


. 


CHAPTER  VII 
CONCLUSIONS 


This  report  has  established  a general  windowing  procedure  for 
discrete-time  feedback  control  system  design.  Initially,  the  basic 
principle  of  Gibbs'  phenomena  was  derived;  and  then  developed  Into 
the  generalized  Gibbs'  phenomena  in  both  the  time  and  frequency  domain. 

In  other  words,  the  reasons  for  using  window  functions  in  compensator 
synthesis  has  been  fully  explored.  The  rationale  were  presented  step 
by  step  In  order  to  develop  a sound  foundation. 

The  process  for  nonrecursive  filter  design  was  presented  and  two 
new  techniques  were  developed.  One  used  the  trapezoidal  windows  to 
offer  a unified  approach,  and  the  other  extended  the  frequency  trans- 
formation method  for  nonrecursive  filter  design. 

In  the  analysis  of  discrete-time  feedback  systems,  we  mainly  applied 
the  Fast  Fourier  Transform  to  avaluate  the  transient  and  frequency 
responses  and  to  study  stabilities.  A new  criterion  based  on  the  return 
difference  was  formulated  which  is  simpler  than  the  Nyquist  criterion. 

Finally,  a detailed  design  procedure  for  discrete-time  feedback 
systems  was  established.  Two  kinds  (type  "0"  and  type  "1")  of  feed- 
back systems  were  used  as  demonstration  examples.  The  five  steps  for 
design  were:  (1)  Specifying  the  required  performance;  (2)  Constructing 

a reference  model;  (3)  Computing  a suitable  frequency  response  of  the 
compensator;  (A»)  Designing  a nonrecursive  filter  transfer  function  by 
the  window  method;  and  (5)  Verifying  the  performance  of  the  compensated 
system.  Compared  with  other  methods,  this  new  window  method  for  compen- 
sator design  is  not  only  simpler  and  more  systematic,  but  also  a fast 
algorithmic  which  Is  computer-oriented. 


APPENDIX  A 

EVALUATION  OF  FOURIER  COEFFICIENTS  AND  COMPUTATION  OF  THE  FOURIER 
TRANSFORM  VIA  THE  FFT  ALGORITHM 

A.O  INTRODUCTION 


The  definitions^  of  the  FFT  algorithm  are 


N_j  -J2irnm 

X(n)  • 2 « N 

mmO 


(A-l) 


and 


N-|  j2wffln 

x(m)  - £ X(n)  e 
n«0 


(A-2) 


where  X(n)  and  x(m)  are  the  sequences  In  the  frequency  domain  and  the 
time  domain,  respectively.  For  the  sake  of  simplicity,  let  us  rewrite 
Eqs.  (A-l)  and  (A-2)  as 


X(n)  - FFMn(x(m)] 


(A-3) 


and 


x(m)  - FF  '(X(n)J 


(A-k) 
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A. 1 EVALUATION  OF  FOURIER  COEFFICIENTS 


A. 1.1  FOURIER  ANALYSIS  OF  PERIODIC  FUNCTIONS  IN  THE  TIME  DOMAIN 


Let  f(t)  be  a periodic  function  In  the  time  domain  with  period  T. 


185 


( 


PAQg 


n d*» 


The  Fourier  series  expansion  is  given  by 


00 

-£ 


f(t)  - } _ e 


jnw.t 


(A-5) 


where 


. fT  -Jnw0t 

*n  “ T J f(t)  e dt 

* A 


(A-6) 


and  u»0  - j . If  we  take  N samples  in  one  period  and  evaluate  Eq.  (A-6) 
by  summation,  we  have 


N-l 


1 . -jnw.mzst 

an  - y / M f (mAt)  e 0 At 


(A-7) 


where  At  ■ ^ . Write  the  finite  partial  sum  of  Eq.  (A-5)  as 


Defining 


and  substituting  Eq.  (A-9)  Into  Eq.  (A-8)  yields 

f(mAt)  - F^(bn] 


(A-8) 


(A-9) 


(A- 10) 
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From  Eqs.  (A-7)  and  (A-8) , we  have 

bn  “ FFNn[f(rU)I  (A'H) 

A. 1.2  FOURIER  ANALYSIS  OF  PERIODIC  FUNCTIONS  IN  THE  FREQUENCY  DOMAIN 


Let  G(w)  be  a periodic  function  In  the  frequency  domain  with  period 
P.  The  Fourier  series  expansion  of  G(w)  Is  given  by 


G(u>)  - X)  An  « 


(A- 12) 


where 


S.-*f 


Jnt0« 

G U)  e dw 


(A- 13) 


In  P 

and  tQ  ■ p—  • Let  Aw  ■ jj-  , and 


n-N 


A N i 

n ■ 0,  j * I 


N u . 

n * y,  ...i  N “ I 


(A- 14) 


By  the  same  reasoning  used  In  subsection  A. 1.1,  we  have 


G(mAu)  - FF^lBjN 


(A- 15) 


end 


FF"’[6(mA«)] 
Bn N 


(A- 16) 
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f 


Let 


F (nAu) 


F - 

n 


. n . 

n ■ 0,  . . . , r-  - 1 


F (n-NAw)  n ■ y,  ....  M - 1 


(A-21) 


f (mAt) 


f - 

m 


w • 0(  • • • f y " 1 


f (HFHAu)  n ■ y,  . . . , N • 1 


(A-22) 


Then,  we  heve 


F - FFU  f T 
n Nn  m 


(A-23) 


, < F„ 

'm T— 


(A-21.) 
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APPENDIX  B 

THE  INVERSE  FOURIER  TRANSFORM  OF  A PERIODIC  FUNCTION  IN  THE 

FREQUENCY  DOMAIN 


Let  C(u)  be  a periodic  function  in  the  frequency  domain  with  per 
iod  u»s . From  Eq.  (3*0  of  Chapter  1 1 , we  can  expand  C(u>)  as  a Fourier 
series  In  the  frequency  domain. 


C(w) 


k“-°° 


A e'JkT“ 

\ e 


where  T - — , and 


(B-i) 


(B-2) 


Therefore,  from  the  definition  of  the  Inverse  Fourier  Transform,  the 
inverse  transform  of  C(u)  Is  given  by 


C(t)  ft  F"1 [C(«)] 


m 

m J2  \ 6(t_kT>  (B-3) 

K"-« 

and  the  Fourier  coefficients  (A^'s)  form  the  time  sequence  of  the  Inverse 
Fourier  Transform  of  C(w). 
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APPENDIX  C 
KEYBOARD  PROGRAMS 


PROGRAM  1 

DEMONSTRATION  OF  THE  GIBBS'  PHENOMENON  IN  THE  TIME  DOMAIN 


1 

L 

0 

A 

BS 

4096 

7 

L 

l 

10 

CL 

0 

13 

X< 

3 

16 

D 

0 

19 

F 

21 

D 

0 

24 

F 

26 

Y 

807 

29 

Y 

80  7 

32 

• 

PROGRAM  2 

DEMONSTRATION  OF  THE  GIBBS'  PHENOMENON  IN  THE  FREQUENCY  DOMAIN 


1 

L 

0 

4 

85 

4096 

7 

L 

1 

10 

CL 

0 

13 

X< 

3 

16 

D 

0 

19 

F 

21 

0 

0 

24 

F 

26 

D 

0 

29 

TP 

0 

32 

D 

0 

35 

Y 

807 

38 

Y 

807 

41 

• 
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PROGRAM  3 

EVALUATION  OF  THE  TIME  RESPONSES  OF  THE  TIME  DOMAIN  TRAPEZOIDAL 
WINDOWS  AND  THE  WINDOW  EFFECTS 


I L 0 


4 

BS 

128 

7 

CL 

0 

10 

CL 

1 

13 

D 

0 

16 

J 

10 

22 

CL 

1 

25 

BS 

4096 

28 

CL 

0 

51 

2048 

33 

t 

0 

64 

37 

F 

0 

40 

: 

0 

64 

44 

0 

2048 

48 

D 

0 

51 

Y 

807 

54 

V 

0 

2036 

2C67 

59 

s 

0 

62 

D 

0 

65 

Y 

807 

68 

w 

0 

2095 

2126 

73 

V 

0 

2180 

2211 

78 

w ' 

0 

1895 

1926 

83 

V 

0 

1960 

1991 

88 

J 

99 

91 

L 

10 

94 

X 

1 

97 

x» 

2 

100 

X 

1 

103 

s 

1 

0 

50 

108 

D 

1 

111 

X 

1 

114 

A- 

2 

117 

* 

0 

-1 

121 

A* 

I 

124 

< 

126 

L 

99 

129 

• 
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PROGRAM  k 


EVALUATION  OF  THE  FREQUENCY  RESPONSE  OF  THE  FREQUENCY  DOMAIN  TRAPEZOIDAL 

WINDOWS  AND  THE  WINDOW  EFFECTS 


1 

L 

0 

4 

BS 

256 

7 

CL 

0 

10 

CL 

1 

13 

D 

0 

16 

J 

10 

19 

D 

0 

22 

CL 

1 

25 

BS 

4096 

28 

CL 

0 

201 

4095 

33 

•» 

0 

100 

37 

r 

0 

40 

V 

0 

0 

3 

45 

D 

0 

48 

D 

0 

0 

500 

53 

Y 

80  7 

56 

D 

0 

59 

X> 

1 

62 

TL 

1 

65 

D 

1 

0 

500 

70 

Y 

807 

73 

D 

1 

0 

160 

78 

J ' 

99 

81 

L 

10 

84 

X 

1 

87 

X> 

2 

90 

X 

1 

93 

s 

1 

0 

100 

98 

D 

1 

101 

X 

1 

104 

A- 

2 

107 

* 

0 

-1 

111 

A+ 

1 

114 

X> 

1 

117 

A- 

50 

120 

* 

0 

-1 

124 

0 

0 

127 

CL 

0 

0 

155 

132 

A* 

1 

135 

«• 

0 

156 

139 

< 

141 

L 

99 

144 

• 

PROGRAM  5 

DESIGN  OF  NONRECURSIVE  FILTERS 


1 

L 

0 

65 

X 

i 

4 

BS 

4096 

68 

X» 

2 

7 

CL 

0 

71 

X 

1 

10 

D 

0 

74 

$ 

1 

0 

13 

F 

0 

79 

D 

1 

16 

t 

0 

4096 

82 

X 

1 

20 

V 

0 

0 

9 

85 

A- 

2 

25 

D 

0 

88 

♦ 

0 

-1 

28 

«• 

0 

4088 

92 

A+ 

1 

32 

D 

0 

95 

X> 

1 

35 

F 

0 

98 

A- 

3 

38 

# 

0 

4096 

101 

♦ 

0 

-1 

42 

TL 

0 

105 

D 

0 

45 

D 

0 

108 

CL 

0 

0 

48 

Y 

80  7 

113 

A+ 

1 

51 

D 

0 

116 

•«* 

0 

56 

54 

Y 

807 

120 

D 

0 

57 

D 

0 

123 

< 

60 

< 

125 

• 

62 

L 

30 

PROGRAM  6 

EVALUATION  OF  THE  FREQUENCY  RESPONSE  OF  DISCRETE-TIME  CONTROL  SYSTEMS 


1 

L 

0 

4 

D 

0 

7 

CL 

0 

10 

CL 

1 

13 

D 

0 

16 

F 

0 

19 

F 

1 

22 

i 

1 

25 

D 

0 

28 

Y 

810 

31 

Y 

807 

34 

Y 

809 

37 

D 

0 

40 

Y 

807 

43 

D 

0 

46 

Y 

810 

49 

Y 

807 

52 

Y 

809 

55 

• 

I 


PROGRAM  7 

EVALUATION  OF  THE  UNIT  IMPULSE  RESPONSE  OF  DISCRETE-TIME  CONTROL  SYSTEMS 


1 L 
4 BS 
7 CL 
10  CL 
13  D 
16  F 
19  F 

22  t 

25  0 
28  F 
31  t 
35  V 
40  . 


0 

1024 

0 

1 

0 

0 

1 

1 

0 

0 

0 1024 

0 0 


63 


PROGRAM  8 

EVALUATION  OF  THE  UNIT  STEP  RESPONSE  OF  DISCRETE-TIME  CONTROL  SYSTEMS 


1 L 
4 BS 
7 CL 
10  CL 
13  D 
16  V 
21  V 
26  F 
29  F 
32  t 
35  D 
38  F 
41  i 
45  V 
50  S 
53  V 
58  . 


0 

1024 

0 

1 


0 

1 

0 

1 

1 

0 


0 

0 


0 

0 1024 

0 0 

0 

0 0 


10 

10 


63 

63 


195 


PROGRAM  9 

CALCULATION  AND  PLOT  OF  THE  PHASE  ANGLE  OF  THE  RETURN  DIFFERENCE 


1 

L 

0 

4 

CL 

0 

7 

CL 

1 

10 

D 

0 

13 

F 

0 

16 

F 

1 

19 

t 

1 

22 

D 

0 

25 

TP 

0 

28 

D 

0 

31 

Y 

807 

34 

Y 

810 

37 

Y 

809 

40 

• 

PROGRAM  10 

FINDING  THE  FREQUENCY  RESPONSE  FROM  THE  UNIT  STEP  RESPONSE 


1 L 

0 

4 BS 

1024 

7 D 

0 

10  V 

0 

0 

63 

15  D 

0 

18  % 

0 

21  V 

0 

0 

63 

26  F 

0 

29  * 

0 

1024 

33  D 

0 

36  Y 

807 

39  Y 

610 

42  Y 

809 

45  D 

0 

48  Y 

807 

51  D 

0 

54  Y 

807 

57  Y 

810 

60  Y 

609 

63  • 

196 


PROGRAM  11 


FINDING  THE  DESIGN  MODEL  AND  ANALYZING  THE  VELOCITY  ERROR  FROM  THE 

SPECIFIED  STEP  RESPONSE 


1 

L 

0 

116 

CL 

0 

0 

4 

BS 

1024 

120 

X» 

9 

7 

CL 

0 

123 

V 

0 

0 

63 

10 

D 

0 

128 

D 

0 

13 

V 

0 
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DESIGN  OF  TYPE  "O"  SYSTEMS 
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PROGRAM  13 

DESIGN  OF  TYPE  "1"  SYSTEMS 
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